1#!/usr/local/bin/python3.8 2# way to run script example 3# python ~/afni/src/python_scripts/align_epi_anat.py \ 4# -anat anat+orig -epi epi_r1+orig -base_epi median -ex_mode dry_run 5# 6# more examples 7 8# align_epi_anat.py -anat anat+orig -epi epi+orig -epi_base 5 9 10# align_epi_anat.py -anat sb23_mpra+orig -epi epi_r03+orig \ 11# -epi2anat -ex_mode dry_run -epi_base 6 -child_epi epi_r??+orig.HEAD \ 12# -anat2epi -epi2anat -tlrc_apar sb23_mpra_at+tlrc -suffix _alx2 13 14import sys 15import copy 16from time import asctime 17 18# AFNI modules 19from afnipy.afni_base import * 20from afnipy.afni_util import * 21from afnipy.option_list import * 22from afnipy.db_mod import * 23from afnipy import ask_me 24 25g_help_string = """ 26 =========================================================================== 27 align_epi_anat.py - align EPI to anatomical datasets or vice versa 28 29 This Python script computes the alignment between two datasets, typically 30 an EPI and an anatomical structural dataset, and applies the resulting 31 transformation to one or the other to bring them into alignment. 32 33 This script computes the transforms needed to align EPI and 34 anatomical datasets using a cost function designed for this purpose. The 35 script combines multiple transformations, thereby minimizing the amount of 36 interpolation applied to the data. 37 38 Basic Usage: 39 align_epi_anat.py -anat anat+orig -epi epi+orig -epi_base 5 40 41 The user must provide EPI and anatomical datasets and specify the EPI 42 sub-brick to use as a base in the alignment. 43 44 Internally, the script always aligns the anatomical to the EPI dataset, 45 and the resulting transformation is saved to a 1D file. 46 As a user option, the inverse of this transformation may be applied to the 47 EPI dataset in order to align it to the anatomical data instead. 48 49 This program generates several kinds of output in the form of datasets 50 and transformation matrices which can be applied to other datasets if 51 needed. Time-series volume registration, oblique data transformations and 52 Talairach (standard template) transformations will be combined as needed 53 and requested (with options to turn on and off each of the steps) in 54 order to create the aligned datasets. 55 56 **Note the intermediate datasets used to compute the alignment are **not** 57 saved unless one of the -save_xxx options is given. This includes 58 skull-stripped, slice timing corrected and volume registered datasets 59 without alignment. These intermediated datasets are normally deleted. 60 See the -save_xxx section below for more information on saving these 61 datasets for future use. 62 63 Depending upon selected options, the script's output contains the following: 64 Datasets: 65 ANAT_al+orig: A version of the anatomy that is aligned to the EPI 66 EPI_al+orig: A version of the EPI dataset aligned to the anatomy 67 EPI_tlrc_al+tlrc: A version of the EPI dataset aligned to a standard 68 template 69 These transformations include slice timing correction and 70 time-series registration by default. 71 72 Transformation matrices: 73 ANAT_al_mat.aff12.1D: matrix to align anatomy to the EPI 74 EPI_al_mat.aff12.1D: matrix to align EPI to anatomy 75 (inverse of above) 76 EPI_vr_al_mat.aff12.1D: matrix to volume register EPI 77 EPI_reg_al_mat.aff12.1D: matrix to volume register and align epi 78 to anatomy (combination of the two 79 previous matrices) 80 EPI_al_tlrc_mat.aff12.1D: matrix to volume register and align epi 81 to anatomy and put into standard space 82 83 Motion parameters from optional volume registration: 84 EPI_tsh_vr_motion.1D: motion parameters from EPI time-series 85 registration (tsh included in name if slice 86 timing correction is also included). 87 88 where the uppercase "ANAT" and "EPI" are replaced by the prefix names 89 of the input datasets, and the suffix can be changed from "_al" as a user 90 option. 91 92 You can use these transformation matrices later to align other datasets: 93 3dAllineate -cubic -1Dmatrix_apply epi_r1_al_mat.aff12.1D \\ 94 -prefix epi_alman epi_r2+orig 95 96 97 The goodness of the alignment should always be assessed visually. 98 Superficially, most of 3dAllineate's cost functions, and those 99 of registration programs from other packages, will produce a plausible 100 alignment based upon a cursory examination but it may not be the best. 101 You need to examine the results carefully if alignment quality is crucial 102 for your analysis. 103 104 In the absence of a gold standard, and given the low contrast of EPI data, 105 it is difficult to judge alignment quality by just looking at the two 106 volumes. This is the case, even when you toggle quickly between one volume 107 and the next, by turning the color overlay off and using the 'u' key in the 108 slice viewer window. To aid with the assessment of alignment, you can use 109 the -AddEdge option or call the @AddEdge script directly. See the help for 110 @AddEdge for more information on that script. 111 112 The default options assume the epi and anat datasets start off fairly close, 113 as is normally the case when the epi dataset closely precedes or follows an 114 anatomical dataset acquisition. If the two data are acquired over separate 115 sessions, or accurate coordinate data is not available in the dataset header 116 (as sometimes occurs for oblique data), various options allow for larger 117 movement including "-cmass cmass", "-big_move","-giant_move", 118 "-ginormous_move", and -align_centers yes". Each of these options 119 is described below. If the datasets do not share the same 120 coordinate space at all, it may be useful to use the "-ginormous_move", 121 "-align_centers" options or run @Align_Centers script first. 122 123 Although this script has been developed primarily for aligning anatomical T1 124 data with EPI BOLD data, it has also been successfully applied for aligning 125 similar modality data together, including T1-SPGR to T1-SPGR, T1-FLAIR 126 to T1-SPGR, EPI to EPI, T1-SPGR at 7T to T1-SPGR at 3T, EPI-rat1 to 127 EPI-rat2, .... If this kind of alignment is required, the default cost 128 function, the Local Pearson Correlation (lpc), is not appropriate. 129 Other cost functions like lpa or nmi have been seen to work well for 130 intra-modality alignment, using the option "-cost lpa". Also see the the 131 dset1 and dset2 options below for functionally equivalent options to the 132 epi and anat options. 133 134 --------------------------------------------- 135 REQUIRED OPTIONS: 136 137 -epi dset : name of EPI dataset 138 -anat dset : name of structural dataset 139 -epi_base : the epi base used in alignment 140 (0/mean/median/max/subbrick#) 141 142 MAJOR OPTIONS: 143 -help : this help message 144 145 -anat2epi : align anatomical to EPI dataset (default) 146 -epi2anat : align EPI to anatomical dataset 147 148 The following options are equivalent to those epi/anat options above 149 except it is assumed the datasets will have similar modalities if 150 either dset1 or dset2 is specified, and the default cost function is 151 changed to 'lpa' instead of 'lpc'. This should reduce confusion when 152 aligning other types of datasets. Most other options that also have 153 names with anat and epi have corresponding dset1 and dset2 options 154 that are exactly equivalent. 155 156 -dset1 : name of dataset1 157 -dset2 : name of dataset2 158 -dset1to2 : align dataset1 to dataset2 159 -dset2to1 : align dataset2 to dataset1 160 161 162 -suffix ssss: append suffix 'sss' to the original anat/epi dataset to use 163 in the resulting dataset names (default is "_al") 164 165 -child_epi dset1 dset2 ... : specify other EPI datasets to align. 166 Time series volume registration will be done to the same 167 base as the main parent EPI dataset. 168 Note if aligning anat to epi, you can still use the -save_vr option 169 to save the volume registered (motion corrected) datasets. See the 170 -save_xxx option section of this help for more information. 171 -child_dset2 equivalent to child_epi above 172 173 -child_anat dset1 dset2 ... : specify other anatomical datasets to align. 174 The same transformation that is computed for the parent anatomical 175 dataset is applied to each of the child datasets. This only makes 176 sense for anat2epi transformations. Skullstripping is not done for 177 the child anatomical dataset. 178 -child_dset1 equivalent to child_anat above 179 180 -AddEdge : run @AddEdge script to create composite edge images of 181 the base epi or anat dataset, the pre-aligned dataset and 182 the aligned dataset. Datasets are placed in a separate 183 directory named AddEdge. The @AddEdge can then be used 184 without options to drive AFNI to show the epi and anat 185 datasets with the edges enhanced. For the -anat2epi case 186 (the default), the anat edges are shown in purple, and the 187 epi edges are shown in cyan (light blue). For the -epi2anat 188 case, the anat edges are shown in cyan, and the epi edges 189 are purple. For both cases, overlapping edges are shown in 190 dark purple. 191 192 -big_move : indicates that large displacement is needed to align the 193 two volumes. This option is off by default. 194 -giant_move : even larger movement required - uses cmass, two passes and 195 very large angles and shifts. May miss finding the solution 196 in the vastness of space, so use with caution 197 -ginormous_move : adds align_centers to giant_move. Useful for very far 198 apart datasets 199 200 Notes on the big_move and giant_move options: 201 "big_move" allows for a two pass alignment in 3dAllineate. 202 The two-pass method is less likely to find a false minimum 203 cost for alignment because it does a number of coarse (blurred, 204 rigid body) alignments first and then follows the best of these 205 coarse alignments to the fine alignment stage. The big_move 206 option should be a relatively safe option, but it adds 207 processing time. 208 209 The giant_move option expands the search parameters in space 210 from 6 degrees and 10 mm to 45 degrees and 45 mm and adds in 211 a center of mass adjustment. The giant_move option will usually 212 work well too, but it adds significant time to the processing 213 and allows for the possibility of a very bad alignment.Another cost 214 functional is available that has worked well with noisy data, "lpc+ZZ". 215 For difficult data, consider that alternative. 216 217 If your data starts out fairly close (probably the typical case 218 for EPI and anatomical data), you can use the -big_move with 219 little problem. All these methods when used with the default 220 lpc cost function require good contrast in the EPI image so that 221 the CSF can be roughly identifiable. 222 223 -rigid_body Limit transformation to translation and rotation, 224 no scaling or shearing. 225 -rigid_equiv Compute alignment with full affine 12 parameters, but 226 use only the translation and rotation parameters. Useful 227 for axialization/AC-PC alignment to a template 228 229 -partial_coverage: indicates that the EPI dataset covers only a part of 230 the brain. Alignment will try to guess which direction should 231 not be shifted If EPI slices are known to be a specific 232 orientation, use one of these other partial_xxxx options. 233 -partial_axial 234 -partial_coronal 235 -partial_sagittal 236 237 -keep_rm_files : keep all temporary files (default is to remove them) 238 -prep_only : do preprocessing steps only 239 -verb nn : provide verbose messages during processing (default is 0) 240 -anat_has_skull yes/no: Anat is assumed to have skull ([yes]/no) 241 -epi_strip methodname : method to mask brain in EPI data 242 ([3dSkullStrip]/3dAutomask/None) 243 -volreg_method methodname: method to do time series volume registration 244 (motion correction) of EPI data 245 ([3dvolreg],3dWarpDrive,3dAllineate). 246 3dvolreg is for 6 parameter (rigid-body) 247 3dWarpDrive is for 12 parameter (general affine) 248 3dAllineate - also 12 parameter with LPA cost function 249 250 Note if aligning anat to epi, the volume registered EPI 251 dataset is **not** saved unless you use the -save_vr 252 option. See the -save_xxx option section of this help for 253 more information. 254 255 -dset1_strip : skull stripping method for dataset1 256 -dset2_strip : skull stripping method for dataset2 (equivalent to epi_strip) 257 258 A template registered anatomical dataset such as a talairach-transformed 259 dataset may be additionally specified so that output data are 260 in template space. The advantage of specifying this transform here is 261 that all transformations are applied simultaneously, thereby minimizing 262 data interpolation. 263 264 -tlrc_apar ANAT+tlrc : structural dataset that has been aligned to 265 a master template such as a tlrc dataset. If this option 266 is supplied, then an epi+tlrc dataset will be created. 267 The @auto_tlrc script may be used to create this 268 "talairach anatomical parent". This option is only valid 269 if aligning epi to anat. 270 271 272 Other options: 273 -ex_mode modename : execute mode (echo/dry_run/quiet/[script]). 274 "dry_run" can be used to show the commands that 275 would be executed without actually running them. 276 "echo" shows the commands as they are executed. 277 "quiet" doesn't display commands at all. 278 "script" is like echo but doesn't show stdout, stderr 279 header lines and "cd" lines. 280 "dry_run" can be used to generate scripts which can be 281 further customized beyond what may be available through 282 the options of this program. 283 -Allineate_opts '-ssss -sss' : options to use with 3dAllineate. Default 284 options are 285 "-weight_frac 1.0 -maxrot 6 -maxshf 10 -VERB -warp aff " 286 -volreg [on]/off : do volume registration on EPI dataset before alignment 287 -volreg_opts '-ssss -sss' : options to use with 3dvolreg 288 -volreg_base basenum/type : the epi base used in time series volume 289 registration. 290 The default is to use the same base as the epi_base. 291 If another subbrick or base type is used, an additional 292 transformation will be computed between the volume 293 registration and the epi_base 294 (0/mean/median/max/subbrick#) 295 296 Note if aligning anat to epi, the volume registered EPI 297 dataset is **not** saved unless you use the -save_vr 298 option. See the -save_xxx option section of this help for 299 more information. 300 301 -tshift [on]/off : do time shifting of EPI dataset before alignment 302 -tshift_opts : options to use with 3dTshift 303 The script will determine if slice timing correction is 304 necessary unless tshift is set to off. 305 306 -deoblique [on]/off : deoblique datasets before alignment 307 -deoblique_opts '-ssss -sss': options to use with 3dWarp deobliquing 308 The script will try to determine if either EPI or anat data 309 is oblique and do the initial transformation to align anat 310 to epi data using the oblique transformation matrices 311 in the dataset headers. 312 313 -master_epi nnn : master grid resolution for aligned epi output 314 -master_tlrc nnn : master grid resolution for epi+tlrc output 315 -master_anat nnn : master grid resolution for aligned anatomical data output 316 -master_dset1 nnn : equivalent to master_anat above 317 -master_dset2 nnn : equivalent to master_epi above 318 (SOURCE/BASE/MIN_DXYZ/dsetname/n.nn) 319 Each of the 'master' options can be set to SOURCE,BASE, 320 a specific master dataset, MIN_DXYZ or a specified cubic 321 voxel size in mm. 322 323 MIN_DXYZ uses the smallest voxel dimension as the basis 324 for cubic output voxel resolution within the bounding box 325 of the BASE dataset. 326 327 SOURCE and BASE are used as in 3dAllineate help. 328 329 The default value for master_epi and master_anat is SOURCE, 330 that is the output resolution and coordinates should be 331 the same as the input. This is appropriate for small 332 movements. 333 334 For cases where either dataset is oblique (and larger 335 rotations can occur), the default becomes MIN_DXYZ. 336 337 The default value for master_tlrc is MIN_DXYZ. 338 339 "-master_dset1" and "-master_dset2" may be used as 340 equivalent expressions for anat and epi output resolutions, 341 respectively. 342 343 -check_flip : check if data may have been left/right flipped by aligning 344 original and flipped versions and then comparing costs 345 between the two. This option produces the L/R flipped 346 and aligned anat/dset1 dataset. A warning is printed 347 if the flipped data has a lower cost function value 348 than the original dataset when both are aligned to the 349 epi/dset2 dataset. 350 351 This issue of left-right confusion can be caused 352 by problems with DICOM files or pipelines 353 that include Analyze format datasets. In these cases, 354 the orientation information is lost, and left-right may 355 be reversed. Other directions can also be confused, but 356 A-P and I-S are usually obvious. Note this problem has 357 appeared on several major publicly available databases. 358 Even if other software packages may proceed without errors 359 despite inconsistent, wrong or even missing coordinate 360 and orientation information, this problem can be easily 361 identified with this option. 362 363 This option does not identify which of the two datasets 364 need to be flipped. It only determines there is likely 365 to be a problem with one or the other of the two input 366 datasets. Importantly, we recommend properly visualizing 367 the datasets in the afni GUI. Look for asymmetries in the 368 two aligned anat/dset1 datasets, and see how they align 369 with the epi/dset2 dataset. To better determine the left 370 and right of each dataset, we recommend relying on tags 371 like vitamin E or looking for surgical markers. 372 373 -flip_giant : apply giant_move options to flipped dataset alignment 374 even if not using that option for original dataset 375 alignment 376 377 -save_xxx options 378 Normally all intermediate datasets are deleted at the end of the script. 379 If aligning anat to epi, the volume registered EPI dataset, although 380 computed, is **not** saved unless you use the -save_vr option. 381 Similarly other intermediate datasets are not saved unless explicitly 382 requested with one of these options: 383 -save_Al_in : save 3dAllineate input files 384 -save_tsh : save tshifted epi 385 -save_vr : save volume registered epi 386 -save_skullstrip : save skull-stripped (not aligned) 387 -save_rep : save representative tstat epi 388 -save_resample : save resampled epi 389 -save_epi_ns : save skull-stripped epi 390 -save_all : save all the above datasets 391 392 Not included with -save_all (since parameters are required): 393 394 -save_orig_skullstrip PREFIX : save original skull-stripped dset 395 -save_script SCRIPT_NAME : save shell command script to given file 396 397 Alternative cost functions and methods: 398 The default method used in this script is the LPC (Localized Pearson 399 Correlation) function. The 'lpc' cost function is computed by the 400 3dAllineate program. Other cost functionals are available and are 401 described briefly in the help for 3dAllineate. This script allows 402 the user to choose any cost function available in that program with 403 404 -cost xxx 405 406 Some cost functionals have proven to be useful for some situations. 407 Briefly, when aligning similar datasets (anat to anat), the 'lpa' method 408 usually provides good alignment. Instead of using a negative correlation, 409 as the 'lpc' method does, the 'lpa' cost functional uses the absolute value 410 of the local correlation, so both positive and negative correlations drive 411 the alignment. Occasionally the simplest least squares cost functional 412 will be useful (implemented with -ls). 413 414 If either of the input datasets has very little structural detail (less 415 than typical EPI), the mutual information methods provide a rough 416 alignment that gives alignment of mostly the contour of the datasets. 417 These are implemented with '-cost nmi' or '-cost mi'. 418 419 The lpa cost function looks for both high positive and negative 420 local Pearson correlation (LPA is an acronym in our program for the 421 absolute value of the local Pearson correlation). The LPC method looks 422 for negative correlation, essentially matching the dark CSF in T1 images 423 with the bright CSF in EPI images. The more negative the correlation the 424 more likely the CSF will overlay each other and carry the rest of the 425 volume along with it. 426 427 -multi_cost cf1 cf2 ... 428 Besides cost from specified cost function or default cost function, 429 also compute alignment using other cost functionals. For example, using 430 "-cost lpa -multi_cost ls nmi" will compute an alignment for the lpa, ls 431 and nmi cost functionals. See 3dAllineate's HELP for a full list of 432 available cost functionals. Use the AFNI GUI to view differences among 433 cost functionals. 434 435 -check_cost cf1 cf2 ... 436 Verify alignment against another cost functional. If there is a large 437 difference, a warning is printed. This does not mean the alignment is 438 bad, only that it is different. 439 440 -edge : use edge method 441 442 The Edge method 443 Finally, the "edge" method is a new method that is implemented not as a 444 cost functional but as a different algorithm altogether. Based on our 445 visualization methods for verifying alignment (as we do in AddEdge), 446 it uses a local approach like the LPA/C cost functionals, but it is 447 independent of the cost function. 448 449 This method has turned out to be useful in a variety of circumstances. It 450 has proven useful for data that changes dramatically over time like 451 manganese-enhanced MRI (MEMRI) and for some data that has other large 452 non-uniformities issues helping to compensate for those large contrasts. 453 454 The edge method prepares the image to be a local spatial variance version 455 of the original image. First both input datasets are automasked with the 456 outer voxel layers removed. The spatial variance is computed over that 457 mask. The optimal alignment is computed between the edge images. Strictly 458 speaking, the datasets are not "edges" but a kind of normalized 2D 459 gradient. The original datasets are then aligned using the transformation 460 computed by the edge image alignment. Internally within the script, 461 the gradient function is accomplished by the 3dLocalstat program using its 462 cvar option for coefficient of variation. The coefficient of variation is 463 computed as the standard deviation within the local neighborhood divided 464 by the mean. The local spatial variance ends up being similar to locally 465 normalized images of edges within the image. 466 467 The "-edge" option is relatively insensitive to most of the cost functions 468 in 3dAllineate, so "lpa", "mi", "nmi" and even "ls" will usually work well. 469 The default is to use the lpa cost functional together with the edge 470 method. 471 472 The edge image is different in a couple ways from the LPA/C correlation. 473 First it is a different function, essentially only a standard deviation 474 over a neighborhood, and then normalized by the absolute value of the 475 mean - effectively a spatial variance (or square root of the variance). 476 The second difference is that while the LPA/C cost functions also operates 477 on local neighborhoods, those neighborhoods are 3-dimensional and set by 478 a neighborhood size set in mm. The shape of the neighborhoods are 479 dodecahedrons (12-side figures) that cover the volume. The edge method 480 instead computes the neighborhoods at each voxel, and the neighborhoods 481 are only two-dimensional - just the voxel and its 8 neighbors in x and y, 482 presumed to be in the same slice rather than across slices. That's for 483 both speed in computation and to remove effects of interpolation or false 484 edges across the relatively distant slices. 485 486 Although not as rigorously tested as the LPC method, this edge method 487 seems to give similar results most of the time. The method does have a few 488 disadvantages compared to the LPC/LPA methods. First, the AddEdge 489 visualization in this script does not support this well (effectively, 490 showing edges of edges). Second, the edge method does not provide 491 three-dimensional edge detection directly. Many times this is an advantage, 492 but if the data has particularly fine slicing in the z-direction, or the 493 data has been resampled, this method may not work as well. Also the method 494 uses an automask to reduce the data so that outside edges do not drive 495 the alignment. The five voxel layer was only empirically found to be 496 useful for this, but may, in fact, be problematic for small partial volumes 497 or for surface coil data where much of the data may be in the area that 498 is masked away. 499 500 The edge method makes no assumption about contrasts between images. Only 501 that edges of features will overlap - the same feature we use visually to 502 verify alignment. This makes it appropriate for both similar and differing 503 modality datasets. 504 505 Both the LPA/LPC and the edge methods require internal features to be 506 present and mostly corresponding in both input datasets. In some cases, 507 this correspondence is not available for aligning some kinds of data with 508 an anatomical references - low-contrast EPI data, radiopharmaceutical PET 509 data targeting specific function, derived parameters from modeling. 510 In these cases, fine alignment is not possible, but alternative cost 511 functions like mutual information or least squares can provide a rough 512 alignment of the contours. 513 514 -output_dir dirname : the default output will put the result in 515 the current directory even if the anat and epi datasets are in other 516 directories. If a directory is specified here, output data including 517 temporary output data will be placed in that directory. If a new directory 518 is specified, that directory will be created first. 519 520 Other obscure and experimental options that should only be handled with 521 care, lest they get out, are visible with -option_help. 522 523 Examples: 524 # align anat to sub-brick 5 of epi+orig. In addition, do slice timing 525 # correction on epi+orig and register all sub-bricks to sub-brick 5 526 # (Sample data files are in AFNI_data4/sb23 in sample class data) 527 # Note the intermediate file, the volume registered EPI dataset, 528 # is **not** saved unless the -save_vr option is also used. 529 # See the -save_xxx option section of this help for more information. 530 531 align_epi_anat.py -anat sb23_mpra+orig -epi epi_r03+orig \\ 532 -epi_base 5 533 534 # Instead of aligning the anatomy to an epi, transform the epi 535 # to match the anatomy. Transform other epi run datasets to be 536 # in alignment with the first epi datasets and with the anatomical 537 # reference dataset. Note that all epi sub-bricks from all runs 538 # are transformed only once in the process, combining volume 539 # registration and alignment to the anatomical dataset in a single 540 # transformation matrix 541 542 align_epi_anat.py -anat sb23_mpra+orig -epi epi_r03+orig \\ 543 -epi_base 5 -child_epi epi_r??+orig.HEAD \\ 544 -epi2anat -suffix al2anat 545 546 # Bells and whistles: 547 # - create Talairach transformed epi datasets (still one transform) 548 # - do not execute, just show the commands that would be executed. 549 # These commands can be saved in a script or modified. 550 # The Talairach transformation requires auto-Talairaching 551 # the anatomical dataset first (cf. @auto_tlrc script) 552 553 @auto_tlrc -base ~/abin/TT_N27+tlrc -input sb23_mpra+orig 554 align_epi_anat.py -anat sb23_mpra+orig -epi epi_r03+orig \\ 555 -epi_base 6 -child_epi epi_r??+orig.HEAD \\ 556 -ex_mode dry_run -epi2anat -suffix _altest \\ 557 -tlrc_apar sb23_mpra_at+tlrc 558 559 560 Our HBM 2008 abstract describing the alignment tools is available here: 561 https://sscc.nimh.nih.gov/sscc/rwcox/abstracts 562 563 Reference: 564 If you find the EPI to Anat alignment capability useful, the paper to 565 cite is: 566 567 ZS Saad, DR Glen, G Chen, MS Beauchamp, R Desai and RW Cox. 568 A new method for improving functional-to-structural alignment using 569 local Pearson correlation. NeuroImage, 44:839-848, 2009. 570 http://dx.doi.org/10.1016/j.neuroimage.2008.09.037 571 572""" 573 574# Also, because input volumes are preprocessed before using 3dAllineate, 575# the script outputs copies of the preprocessed volumes as they were used 576# in 3dAllineate. 577# 578# EPI_epi_in_3dAl_al+orig : EPI volume for 3dAllineate's -base 579# ANAT_epi_in_3dAl_al+orig: ANAT volume for 3dAllineate's -input 580# EPI_wt_3dAl_al+orig : weight volume for 3dAllineate's -weight 581 582# -cost : cost function used by 3dAllineate. Default is lpc, anything 583# else is inferior! 584# -cmass cmass+ss: center of mass option for 3dAllineate 585# ('cmass+a','cmass+xy','nocmass',...) Default is cmass+a. 586# -child_anat dset1 dset2 ... : specify other anatomical datasets to align. 587# The anatomical data will be aligned first to the parent 588# structural dataset. If aligning to EPI data, then the 589# transformation of the parent anatomical to the EPI data will 590# be combined with the inter-structural transformation. 591# 592# -fresh : remove any temporary files at start from a previous run 593# Weighting mask options: 594# A weighting mask is used in the alignment. The default weighting mask 595# is the stripped epi dataset normalized from 0 to 1. 596# -pow_mask n.n : raise the epi masked dataset to a power before normalizing 597# (default is 1.0) 598 599# -tlrc_epar epi_template_dset : EPI dataset that has been aligned to 600# a master template such as a tlrc dataset If this option 601# is supplied, then an anat+tlrc dataset will be created. 602 603 604 605## BEGIN common functions across scripts (loosely of course) 606class RegWrap: 607 def __init__(self, label): 608 self.align_version = "1.62" # software version (update for changes) 609 self.label = label 610 self.valid_opts = None 611 self.user_opts = None 612 self.verb = 1 # a little talkative by default 613 self.save_script = '' # save completed script into given file 614 self.rewrite = 0 #Do not recreate existing volumes 615 self.oexec = "" #dry_run is an option 616 self.epi2anat = 0 # align epi to anat optionally 617 self.anat2epi = 1 # align anat to epi by default 618 self.rmrm = 1 # remove temporary files 619 self.prep_only = 0 # do preprocessing only 620 self.odir = os.getcwd() 621 self.align_centers = 0 # don't align centers 622 self.tshift_flag = 1 # do time shifting on EPI 623 self.volreg_flag = 0 # don't do volume registration on EPI by default 624 self.resample_flag = 1 # do resample 625 self.deoblique_flag = 1 # deoblique datasets first 626 self.deoblique_opt = "" # deobliquing/obliquing options 627 self.skullstrip_opt = "" # skullstripping options 628 self.epistrip_opt = "" # similar stripping/automask options for EPI/dset2 629 self.cmass = "nocmass" # no center of mass option for 3dAllineate 630 self.epi_base = None # don't assume representative epi 631 self.reg_mat = "" # volume registration matrix 1D file 632 self.obl_a2e_mat = "" # oblique anat to epi matrix 633 self.edge = 0 # don't use edges for alignment 634 self.edgelevels = 5 # number of outer layers to remove in edge method 635 self.cost = '' # assign cost below 636 self.epi_dir = '' # assign path from EPI (dset2) dataset's path 637 self.anat_dir = '' # assign path from anat (dset1) dataset's path 638 self.output_dir = '' # user assigned path for anat and EPI 639 self.flip = 0 # don't test for left/right flipping 640 self.flip_giant = 0 # don't necessarily use giant_move for L-R flipping test 641 self.giant_move = 0 # don't use giant_move 642 self.supersize = 0 # don't use supersize 643 self.rigid_equiv = 0 # don't do rigid_body equivalent for alignment 644 645 # options for saving temporary datasets permanently 646 self.save_Al_in = 0 # don't save 3dAllineate input files 647 self.save_tsh = 0 # don't save tshifted epi 648 self.save_vr = 0 # don't save volume registered epi 649 self.save_skullstrip = 0 # don't save skullstripped (not aligned) 650 self.save_origstrip = '' # prefix for simple skullstripped dset 651 self.save_rep = 0 # don't save representative tstat epi 652 self.save_resample = 0 # don't save resampled epi 653 self.save_epi_ns = 0 # don't save skull-stripped epi 654 655 self.master_epi_option = '-master SOURCE' 656 self.master_tlrc_option = '-master SOURCE' # changed to min dimension below 657 self.master_anat_option = '' 658 659 self.dset1_generic_name = 'anat' 660 self.dset2_generic_name = 'epi' 661 662 return 663 664# box, bin and fat mask are not used for now 665 666 def init_opts(self): 667 668 self.valid_opts = OptionList('init_opts') 669 670 self.valid_opts.add_opt('-epi', 1, [], \ 671 helpstr="EPI dataset to align or to which to align") 672 self.valid_opts.add_opt('-dset2', 1, [], \ 673 helpstr="dataset to align or to which to align") 674 675 self.valid_opts.add_opt('-anat', 1, [], \ 676 helpstr="Anatomical dataset to align or to which to align") 677 self.valid_opts.add_opt('-dset1', 1, [], \ 678 helpstr="Dataset to align or to which to align") 679 680 self.valid_opts.add_opt('-keep_rm_files', 0, [], \ 681 helpstr="Don't delete any of the temporary files created here") 682 self.valid_opts.add_opt('-prep_only', 0, [], \ 683 helpstr="Do preprocessing steps only without alignment") 684 self.valid_opts.add_opt('-help', 0, [], \ 685 helpstr="The main help describing this program with options") 686 self.valid_opts.add_opt('-limited_help', 0, [], \ 687 helpstr="The main help without all available options") 688 self.valid_opts.add_opt('-option_help', 0, [], \ 689 helpstr="Help for all available options") 690 self.valid_opts.add_opt('-version', 0, [], \ 691 helpstr="Show version number and exit") 692 self.valid_opts.add_opt('-ver', 0, [], \ 693 helpstr="Show version number and exit") 694 self.valid_opts.add_opt('-verb', 1, [], \ 695 helpstr="Be verbose in messages and options" ) 696 # 26 Nov 2012 [rickr] 697 self.valid_opts.add_opt('-save_script', 1, [], \ 698 helpstr="save executed script in given file" ) 699 700 self.valid_opts.add_opt('-align_centers', 1, ['no'], ['yes', 'no', 'on', 'off'], \ 701 helpstr="align centers of datasets based on spatial\n" \ 702 "extents of the original volume") 703 self.valid_opts.add_opt('-anat_has_skull', 1, [], ['yes', 'no'],\ 704 helpstr="Do not skullstrip anat dataset") 705 self.valid_opts.add_opt('-epi_strip', 1, [], \ 706 ['3dSkullStrip', '3dAutomask', 'None'], \ 707 helpstr="Method to remove skull for EPI data") 708 self.valid_opts.add_opt('-dset1_strip', 1, [], \ 709 ['3dSkullStrip', '3dAutomask', 'None'], \ 710 helpstr="Method to remove skull for dset1 data") 711 self.valid_opts.add_opt('-dset2_strip', 1, [], \ 712 ['3dSkullStrip', '3dAutomask', 'None'], \ 713 helpstr="Method to remove skull for dset2 data") 714 715 self.valid_opts.add_opt('-volreg_method', 1, ['3dvolreg'], \ 716 ['3dvolreg', '3dWarpDrive', '3dAllineate'], \ 717 helpstr="Time series volume registration method\n" \ 718 "3dvolreg: rigid body least squares\n" \ 719 "3dWarpDrive: 12 parameter least squares\n" \ 720 "3dAllineate: 12 parameter LPA cost function\n") 721 self.valid_opts.add_opt('-ex_mode', 1, ['script'], \ 722 ['quiet', 'echo', 'dry_run', 'script'], \ 723 helpstr="Command execution mode.\n" \ 724 "quiet: execute commands quietly\n" \ 725 "echo: echo commands executed\n" \ 726 "dry_run: only echo commands\n" ) 727 728 self.valid_opts.add_opt('-overwrite', 0, [],\ 729 helpstr="Overwrite existing files") 730 self.valid_opts.add_opt('-big_move', 0, [], \ 731 helpstr="Large movement between epi and anat.\n" \ 732 "Uses twopass option for 3dAllineate.\n" \ 733 "Consider cmass options, giant_move,\n" \ 734 "ginormous_move or -align_centers") 735 self.valid_opts.add_opt('-giant_move', 0, [], \ 736 helpstr="Even larger movement between epi and anat.\n" \ 737 "Uses twopass option for 3dAllineate.\n" \ 738 "cmass options and wide angles and shifts") 739 self.valid_opts.add_opt('-ginormous_move', 0, [], \ 740 helpstr="Adds align_centers to giant_move") 741 self.valid_opts.add_opt('-supersize', 0, [], \ 742 helpstr="Large scaling difference - up to 50%") 743 744 745 self.valid_opts.add_opt('-rigid_body', 0, [], \ 746 helpstr="Do only rigid body alignment - shifts and rotates") 747 self.valid_opts.add_opt('-rigid_equiv', 0, [], \ 748 helpstr="Do only rigid body equivalent alignment - shifts and rotates") 749 750 self.valid_opts.add_opt('-partial_coverage', 0, [], \ 751 helpstr="partial_xxxx options control center of mass adjustment") 752 self.valid_opts.add_opt('-partial_axial', 0, []) 753 self.valid_opts.add_opt('-partial_coronal', 0, []) 754 self.valid_opts.add_opt('-partial_sagittal', 0, []) 755 self.valid_opts.add_opt('-AddEdge', 0, [], helpstr= \ 756 "Run @AddEdge script to create double-edge images") 757 758 self.valid_opts.add_opt('-Allineate_opts', -1, \ 759 ["-weight_frac 1.0 -maxrot 6 -maxshf 10 -VERB"\ 760 " -warp aff -source_automask+4 "],\ 761 helpstr="Options passed to 3dAllineate.") 762 763 764 self.valid_opts.add_opt('-perc', 1, ['90']) 765# self.valid_opts.add_opt('-fresh', 0, []) 766 self.valid_opts.add_opt('-suffix', 1,['_al']) 767 self.valid_opts.add_opt('-cost', 1,[]) 768# self.valid_opts.add_opt('-fat', 1, ['1']) 769 self.valid_opts.add_opt('-multi_cost', -1,[], \ 770 helpstr = "can use multiple cost functionals (lpc,lpa,nmi,....\n" \ 771 "See 3dAllineate -HELP for the full list\n") 772 self.valid_opts.add_opt('-check_cost', -1,[], \ 773 helpstr = "Verify alignment against another method\n" 774 "Can use multiple cost functionals (lpc,lpa,nmi,....\n" \ 775 "See 3dAllineate -HELP for the full list\n") 776 777 # transform anat to epi by default, but allow the other way 778 # the resulting transformation will be done at the end to include 779 # any volreg and oblique transformations 780 self.valid_opts.add_opt('-epi2anat', 0, [], \ 781 helpstr = "align EPI dataset to anat dataset") 782 self.valid_opts.add_opt('-anat2epi', 0, [], \ 783 helpstr = "align anat dataset to EPI dataset (default)") 784 self.valid_opts.add_opt('-dset2to1', 0, [], \ 785 helpstr = "align dset2 dataset to dset1 dataset") 786 self.valid_opts.add_opt('-dset1to2', 0, [], \ 787 helpstr = "align dset1 dataset to dset2 dataset (default)") 788 789 # select base EPI dataset type 790 self.valid_opts.add_opt('-epi_base', 1, [], [], \ 791 helpstr = "Base sub-brick to use for alignment\n" \ 792 "Choose sub-brick number or statistic type\n" \ 793 "Valid choices can be, for example, 0,5,mean") 794# ['0', 'sub-brick-n', 'mean', 'median', 'max']) 795 self.valid_opts.add_opt('-dset2_base', 1, [], [], \ 796 helpstr = "Base sub-brick to use for alignment\n" \ 797 "Choose sub-brick number or statistic type\n" \ 798 "Valid choices can be, for example, 0,5,mean") 799 800 # select base EPI type for volume registration 801 self.valid_opts.add_opt('-volreg_base', 1, [], [], \ 802 helpstr = "Base to use for volume registration\n" \ 803 "Choose sub-brick number or statistic type\n" \ 804 "Valid choices can be, for example, 0,5,median") 805 806 # do volume registration of EPI as part of this whole mess 807 self.valid_opts.add_opt('-volreg', 1, [], ['on','off']) 808 self.valid_opts.add_opt('-volreg_opts', -1, ["-cubic"]) 809 810 # do time shifting 811 self.valid_opts.add_opt('-tshift', 1, [], ['on','off']) 812 self.valid_opts.add_opt('-tshift_opts', -1, []) 813 814 # obliquity options 815 self.valid_opts.add_opt('-deoblique', 1, [], ['on','off']) 816 self.valid_opts.add_opt('-deoblique_opts', -1, []) 817 818 # resampling epi to anat 819 self.valid_opts.add_opt('-resample', 1, [], ['on', 'off']) 820 821 # turn off all pre-processing steps 822 self.valid_opts.add_opt('-prep_off', 0, [], \ 823 helpstr = "turn off all pre-processing steps including\n" \ 824 "deoblique, tshift, volreg and resample") 825 # 3dAllineate cmass options 826 self.valid_opts.add_opt('-cmass', 1, [], [], \ 827 helpstr = "choose center of mass options for 3dAllineate\n" 828 "Center of mass shifts the center of the datasets to match\n" 829 "by computing the weighted centers of each.\n" 830 "For partial data, this may be too far in one direction\n" 831 "See 3dAllineate help for details\n" 832 "Valid options include cmass+a, cmass+xy, nocmass\n" 833 "nocmass = no center of mass shift - default\n" 834 "cmass = center of mass shift - used with giant, ginormous_move\n" 835 "cmass+a = automatic center of mass for partial data\n" 836 "cmass+xy,xz,yz = automatic center of mass for partial\n" 837 " axial,coronal,sagittal\n" 838 "For partial data, it may be easier to select one\n" 839 " of the partial_... options above" ) 840 841 # talairach transformed anatomical parent dataset 842 self.valid_opts.add_opt('-tlrc_apar', 1, [], \ 843 helpstr="If this is set, the results will include +tlrc\n" 844 "template transformed datasets for the epi aligned\n" 845 "to the anatomical combined with this additional\n" 846 "transformation to template of this parent dataset\n" 847 "The result will be EPI_al+tlrc.HEAD\n") 848 849 # talairach transformed EPI parent dataset 850 self.valid_opts.add_opt('-tlrc_epar', 1, [], \ 851 helpstr="Not available yet.\n" 852 "If this is set, the results will include +tlrc\n" 853 "template transformed datasets for the anatomical\n" 854 "aligned to the epi combined with this additional\n" 855 "transformation to template of this parent dataset\n" 856 "The result will be ANAT_al+tlrc.HEAD\n") 857 858 # auto_talairach results 859 self.valid_opts.add_opt('-auto_tlrc', 0, [], \ 860 helpstr="Not available yet.\n" 861 "If this is set, the results will also be aligned\n" 862 "to a template using the @auto_tlrc script.\n" 863 "Transformations computed from that will be combined\n" 864 "with the anat to epi transformations and epi to anat\n" 865 "(and volreg) transformations\n" 866 "Only one of the -tlrc_apar, -tlrc_epar or the \n" 867 "-auto_tlrc options may be used\n") 868 # child epi datasets 869 self.valid_opts.add_opt('-child_epi', -1,[],\ 870 helpstr="Names of child EPI datasets") 871 self.valid_opts.add_opt('-child_dset2', -1,[],\ 872 helpstr="Names of children of dset2 datasets") 873 874 # child anat datasets 875 self.valid_opts.add_opt('-child_anat', -1,[],\ 876 helpstr="Names of child anatomical datasets") 877 self.valid_opts.add_opt('-child_dset1', -1,[],\ 878 helpstr="Names of children of dset1 datasets") 879 880 # master resampling options for alignment 881 self.valid_opts.add_opt('-master_epi', 1,[],\ 882 helpstr="-master grid resolution for epi to anat alignment\n" 883 "MIN_DXYZ uses the smallest dimension\n" 884 "Other options are SOURCE and BASE as in 3dAllineate\n" 885 "help. For cases where either dataset is oblique, the\n" 886 "default becomes MIN_DXYZ") 887 self.valid_opts.add_opt('-master_dset2', 1,[],\ 888 helpstr="-master grid resolution for epi to anat alignment\n" 889 "MIN_DXYZ uses the smallest dimension\n" 890 "Other options are SOURCE and BASE as in 3dAllineate\n" 891 "help. For cases where either dataset is oblique, the\n" 892 "default becomes MIN_DXYZ") 893 894 self.valid_opts.add_opt('-master_tlrc', 1,[],\ 895 helpstr="-master grid resolution for epi to tlrc anat\n" 896 "alignment\n" 897 "MIN_DXYZ uses the smallest dimension\n" 898 "Other options are SOURCE and BASE as in 3dAllineate\n" 899 "help") 900 901 self.valid_opts.add_opt('-master_anat', 1,[],\ 902 helpstr="-master grid resolution for anat to epi output\n" 903 "MIN_DXYZ uses the smallest dimension\n" 904 "Other options are SOURCE, BASE, 'n' mm or gridset") 905 self.valid_opts.add_opt('-master_dset1', 1,[],\ 906 helpstr="-master grid resolution for dset1 to dset2 output\n" 907 "MIN_DXYZ uses the smallest dimension\n" 908 "Other options are SOURCE, BASE, 'n' mm or gridset") 909 910 self.valid_opts.add_opt('-master_anat_dxyz', -1,[],\ 911 helpstr="-master grid resolution size (cubic only)\n") 912 self.valid_opts.add_opt('-master_dset1_dxyz', -1,[],\ 913 helpstr="-master grid resolution size (cubic only)\n") 914 915 self.valid_opts.add_opt('-master_epi_dxyz', -1,[],\ 916 helpstr="-master grid resolution (cubic only)\n") 917 self.valid_opts.add_opt('-master_dset2_dxyz', -1,[],\ 918 helpstr="-master grid resolution (cubic only)\n") 919 920 self.valid_opts.add_opt('-master_tlrc_dxyz', -1,[],\ 921 helpstr="-master grid resolution (cubic only)\n") 922 923 924 # apply pre/post-transformation matrice 925 self.valid_opts.add_opt('-pre_matrix', 1, [], \ 926 helpstr="Apply an initial transformation from a 1D file (NB:\n" 927 "not from a *.aff12.1D file); the *.1D file should\n" 928 "contain a 3x4 matrix of numbers.\n" 929 "For example, this file may be one generated by \n" 930 "@Align_Centers, or if inverting a matrix, with:\n" 931 " cat_matvec mat.aff12.1D -I > mat_INV.1D\n" 932 "The transformation will be applied to the\n" 933 "anatomical data before aligning to the EPI\n" 934 "instead of using the built-in obliquity matrices,\n" 935 "if any") 936 self.valid_opts.add_opt('-post_matrix', 1, [], \ 937 helpstr="Apply an additional transformation from a 1D file.\n" 938 "This transformation will be applied to the anatomical\n" 939 "data after alignment with the EPI. This will be\n" 940 "applied similarly to the tlrc transformation and in\n" 941 "place of it.\n" 942 "Output datasets are kept in the 'orig' view") 943 self.valid_opts.add_opt('-skullstrip_opts', -1, [], \ 944 helpstr="Alternate options for 3dSkullstrip.\n" 945 "like -rat or -blur_fwhm 2") 946 self.valid_opts.add_opt('-dset1strip_opts', -1, [], \ 947 helpstr="Alternate name for skullstrip_opts") 948 self.valid_opts.add_opt('-epistrip_opts', -1, [], \ 949 helpstr="Alternate options for 3dSkullstrip/3dAutomask.\n" 950 "like -rat or -blur_fwhm 2 or -peels 2") 951 self.valid_opts.add_opt('-dset2strip_opts', -1, [], \ 952 helpstr="Alternate name for epistrip_opts") 953 self.valid_opts.add_opt('-feature_size', 1, [],\ 954 helpstr="Minimal size in mm of structures in images to match.\n"\ 955 "Changes options for 3dAllineate for the coarse\n" \ 956 "blurring and lpc/lpa neighborhood sizes.May be useful\n" \ 957 "for rat brains, anat to anat and other\n" \ 958 "'non-standard' alignment") 959 self.valid_opts.add_opt('-rat_align', 0, [],\ 960 helpstr="Set options appropriate for rat data - \n" 961 "namely skullstrip and feature size options above.\n") 962 self.valid_opts.add_opt('-output_dir', 1,[],\ 963 helpstr="Set directory for output datasets\n") 964 965 # create edge images 966 # do edge-based alignment 967 self.valid_opts.add_opt('-edge', 0, [],\ 968 helpstr="Use internal edges to do alignment") 969 self.valid_opts.add_opt('-edge_erodelevel', 1, [],\ 970 helpstr="Number of layers to remove for edge method") 971 972 self.valid_opts.trailers = 0 # do not allow unknown options 973 974 self.valid_opts.add_opt('-check_flip', 0, [], \ 975 helpstr="Check if L/R flipping gives better results") 976 self.valid_opts.add_opt('-flip_giant', 0, [], \ 977 helpstr="use giant_move on flipped data even if not used\n"\ 978 "on original data") 979 980 # saving optional output datasets 981 # save datasets used as input to 3dAllineate 982 self.valid_opts.add_opt('-save_Al_in', 0, [], \ 983 helpstr = "Save datasets used as input to 3dAllineate") 984 # save vr dataset 985 self.valid_opts.add_opt('-save_vr', 0, [], \ 986 helpstr = "Save motion-corrected epi dataset") 987 # save timeshifted dataset 988 self.valid_opts.add_opt('-save_tsh', 0, [], \ 989 helpstr = "Save time-series corrected dataset") 990 # save skullstripped anat before alignment 991 self.valid_opts.add_opt('-save_skullstrip', 0, [], \ 992 helpstr = "Save unaligned, skullstripped dataset") 993 # save skullstripped anat before alignment and without oblique xform 994 self.valid_opts.add_opt('-save_orig_skullstrip', 1, [], \ 995 helpstr = "Save simply skullstripped dataset") 996 # save skullstripped epi before alignment 997 self.valid_opts.add_opt('-save_epi_ns', 0, [], \ 998 helpstr = "Save unaligned, skullstripped EPI dataset") 999 # save representative epi tstat epi before alignment 1000 self.valid_opts.add_opt('-save_rep', 0, [], \ 1001 helpstr = "Save unaligned representative tstat EPI dataset") 1002 # save resampled epi dataset before alignment 1003 self.valid_opts.add_opt('-save_resample', 0, [], \ 1004 helpstr = "Save unaligned EPI dataset resampled to anat grid") 1005 1006 # save all the optional datasets 1007 self.valid_opts.add_opt('-save_all', 0, [], \ 1008 helpstr = "Save all optional datasets") 1009 1010 1011 # weighting mask options 1012 self.valid_opts.add_opt('-pow_mask', 1, ['1.0'], \ 1013 helpstr = "power for weighting 1 or 2") 1014 self.valid_opts.add_opt('-bin_mask', 1, ['no'], ['yes', 'no'], \ 1015 helpstr = "convert weighting mask to 0 or 1 - Unused") 1016 self.valid_opts.add_opt('-box_mask', 1, ['no'], ['yes', 'no'], \ 1017 helpstr = "Unused") 1018 1019 self.valid_opts.add_opt('-mask', -1, ['vent'], \ 1020 helpstr="Not available yet.\n" 1021 "Mask to apply to data.") 1022 1023 1024 1025 def dry_run(self): 1026 if self.oexec != "dry_run": 1027 return 0 1028 else: 1029 return 1 1030 1031 def apply_initial_opts(self, opt_list): 1032 opt1 = opt_list.find_opt('-version') # user only wants version 1033 opt2 = opt_list.find_opt('-ver') 1034 if ((opt1 != None) or (opt2 != None)): 1035 # ps.version() 1036 ps.ciao(0) # terminate 1037 opt = opt_list.find_opt('-verb') # set and use verb 1038 if opt != None: self.verb = int(opt.parlist[0]) 1039 1040 opt = opt_list.find_opt('-save_script') # save executed script 1041 if opt != None: self.save_script = opt.parlist[0] 1042 1043 opt = opt_list.find_opt('-ex_mode') # set execute mode 1044 if opt != None: self.oexec = opt.parlist[0] 1045 1046 opt = opt_list.find_opt('-keep_rm_files') # keep temp files 1047 if opt != None: self.rmrm = 0 1048 1049 opt = opt_list.find_opt('-prep_only') # preprocessing only 1050 if opt != None: self.prep_only = 1 1051 1052 opt = opt_list.find_opt('-help') # does the user want help? 1053 if opt != None: 1054 ps.self_help(2) # always give full help now by default 1055 ps.ciao(0) # terminate 1056 1057 opt = opt_list.find_opt('-limited_help') # less help? 1058 if opt != None: 1059 ps.self_help() 1060 ps.ciao(0) # terminate 1061 1062 opt = opt_list.find_opt('-option_help') # help for options only 1063 if opt != None: 1064 ps.self_help(1) 1065 ps.ciao(0) # terminate 1066 1067 opt = opt_list.find_opt('-perc') # set and use percentile for weight 1068 if opt != None: self.perc = float(opt.parlist[0]) 1069 1070 opt = opt_list.find_opt('-suffix') 1071 if opt != None: 1072 self.suffix = opt.parlist[0] 1073 if((opt=="") or (opt==" ")) : 1074 self.error_msg("Cannot have blank suffix") 1075 ps.ciao(1); 1076 1077 opt = opt_list.find_opt('-cost') 1078 if opt != None: self.cost = opt.parlist[0] 1079 else: self.cost = '' 1080 1081 opt = opt_list.find_opt('-pow_mask') 1082 if opt != None: self.sqmask = opt.parlist[0] 1083 1084 opt = opt_list.find_opt('-box_mask') 1085 if opt != None: self.boxmask = opt.parlist[0] 1086 1087 opt = opt_list.find_opt('-bin_mask') 1088 if opt != None: self.binmask = opt.parlist[0] 1089 1090 if (opt_list.find_opt('-epi2anat') or 1091 opt_list.find_opt('-dset2to1') ): # align epi to anat 1092 self.epi2anat = 1 1093 self.anat2epi = 0 # turn off anat to epi unless requested 1094 self.volreg_flag = 1 # turn on motion correction/volume registration 1095 if(opt_list.find_opt('-anat2epi') or 1096 opt_list.find_opt('-dset1to2')): # align anat to epi 1097 self.anat2epi = 1 1098 1099 opt = opt_list.find_opt('-prep_off') # turn off all preprocessing steps 1100 if opt != None: 1101 self.deoblique_flag = 0 1102 self.tshift_flag = 0 1103 self.volreg_flag = 0 1104 self.resample_flag = 0 1105 self.info_msg( \ 1106 "turning off deobliquing tshift, volume registration, resampling") 1107 # note - individual flags can be turned on, if desired 1108 1109 opt = opt_list.find_opt('-deoblique') # deoblique data 1110 if opt != None: 1111 if(opt.parlist[0]=='on'): 1112 self.deoblique_flag = 1 1113 elif(opt.parlist[0]=='off'): 1114 self.deoblique_flag = 0 1115 else: 1116 self.error_msg("deoblique option not on/off") 1117 self.ciao(1) 1118 1119 opt = opt_list.find_opt('-tshift') # do time shifting 1120 if opt != None: 1121 if(opt.parlist[0]=='on'): 1122 self.tshift_flag = 1 1123 elif(opt.parlist[0]=='off'): 1124 self.tshift_flag = 0 1125 else: 1126 self.error_msg("tshift option not on/off") 1127 self.ciao(1) 1128 1129 opt = opt_list.find_opt('-volreg') # do volume registration 1130 if opt != None: 1131 if(opt.parlist[0]=='on'): 1132 self.volreg_flag = 1 1133 self.info_msg("turning on volume registration") 1134 elif(opt.parlist[0]=='off'): 1135 self.volreg_flag = 0 1136 self.info_msg("turning off volume registration") 1137 else: 1138 self.error_msg("volreg option not on/off"); 1139 self.ciao(1) 1140 1141 opt = opt_list.find_opt('-resample') # resample epi to anat 1142 if opt != None: 1143 if(opt.parlist[0]=='on'): 1144 self.resample_flag = 1 1145 elif(opt.parlist[0]=='off'): 1146 self.resample_flag = 0 1147 self.info_msg("turning off resampling") 1148 else: 1149 self.error_msg("resample option not on/off"); 1150 self.ciao(1) 1151 1152 opt = opt_list.find_opt('-check_flip') # check for left/right flipping 1153 if opt != None: 1154 self.flip = 1 1155 self.save_epi_ns = 1 # save EPI for QC, even if anat to epi 1156 1157 # optional data to save 1158 opt = opt_list.find_opt('-save_Al_in') # save 3dAllineate input datasets 1159 if opt != None: self.save_Al_in = 1 1160 opt = opt_list.find_opt('-save_tsh') # save timeshifted data 1161 if opt != None: self.save_tsh = 1 1162 opt = opt_list.find_opt('-save_vr') # save volume registered epi data 1163 if opt != None: 1164 self.save_vr = 1 1165 opt = opt_list.find_opt('-volreg') # if volreg has not been set 1166 if opt == None: self.volreg_flag = 1 # turn on volreg processing 1167 opt = opt_list.find_opt('-save_skullstrip') # save unaligned skullstripped 1168 if opt != None: self.save_skullstrip = 1 1169 # save unaligned, unobliqued dataset 1170 opt = opt_list.find_opt('-save_orig_skullstrip') 1171 if opt != None: 1172 val, err = opt_list.get_string_opt('', opt=opt) 1173 if val == None or err: 1174 ps.self_help() 1175 ps.ciao(0) # terminate 1176 self.save_origstrip = val 1177 opt = opt_list.find_opt('-save_epi_ns') # save unaligned skullstripped epi 1178 if opt != None: self.save_epi_ns = 1 1179 opt = opt_list.find_opt('-save_rep') # save unaligned representative epi 1180 if opt != None: self.save_rep = 1 1181 opt = opt_list.find_opt('-save_resample') # save unaligned resampled epi 1182 if opt != None: self.save_resample = 1 1183 1184 opt = opt_list.find_opt('-save_all') # save all optional output datasets 1185 if ((opt != None) or self.prep_only): 1186 self.save_Al_in = 1 # save 3dAllineate input files 1187 self.save_tsh = 1 # save tshifted epi 1188 self.save_vr = 1 # save volume registered epi 1189 self.save_skullstrip = 1 # save skullstripped (not aligned) 1190 self.save_rep = 1 # save representative tstat epi 1191 self.save_resample = 1 # save resampled epi 1192 self.save_epi_ns = 1 # save skull-stripped epi 1193 1194 1195 1196 def get_user_opts(self): 1197 self.valid_opts.check_special_opts(sys.argv) #ZSS March 2014 1198 self.user_opts = read_options(sys.argv, self.valid_opts) 1199 if self.user_opts == None: return 1 #bad 1200 # no options: apply -help 1201 if ( len(self.user_opts.olist) == 0 or \ 1202 len(sys.argv) <= 1 ) : 1203 ps.self_help() 1204 ps.ciao(0) # terminate 1205 if self.user_opts.trailers: 1206 opt = self.user_opts.find_opt('trailers') 1207 if not opt: 1208 print("** ERROR: seem to have trailers, but cannot find them!") 1209 else: 1210 print("** ERROR: have invalid trailing args: %s", opt.show()) 1211 return 1 # failure 1212 1213 # apply the user options 1214 if self.apply_initial_opts(self.user_opts): return 1 1215 1216 if self.verb > 3: 1217 self.show('------ found options ------ ') 1218 1219 return 1220 1221 def show(self, mesg=""): 1222 print('%s: %s' % (mesg, self.label)) 1223 if self.verb > 2: self.valid_opts.show('valid_opts: ') 1224 self.user_opts.show('user_opts: ') 1225 1226 def info_msg(self, mesg=""): 1227 if(self.verb >= 1) : 1228 print("#++ %s" % mesg) 1229 1230 def error_msg(self, mesg=""): 1231 print("#**ERROR %s" % mesg) 1232 1233 def exists_msg(self, dsetname=""): 1234 print("** Dataset: %s already exists" % dsetname) 1235 print("** Not overwriting.") 1236 if(not ps.dry_run()): 1237 self.ciao(1) 1238 1239 def ciao(self, i): 1240 if i > 0: 1241 print("** ERROR - script failed") 1242 elif i==0: 1243 print("") 1244 1245 os.chdir(self.odir) 1246 1247 if self.save_script: 1248 write_afni_com_history(self.save_script) 1249 1250 # return status code 1251 sys.exit(i) 1252 1253 # save the script command arguments to the dataset history 1254 def save_history(self, dset, exec_mode): 1255 self.info_msg("Saving history") # sounds dramatic, doesn't it? 1256 cmdline = args_as_command(sys.argv, \ 1257 '3dNotes -h "', '" %s' % dset.input()) 1258 com = shell_com( "%s\n" % cmdline, exec_mode) 1259 com.run() 1260 1261 # show help 1262 # if help_level is 1, then show options help only 1263 # if help_level is 2, then show main help and options help 1264 def self_help(self, help_level=0): 1265 if(help_level!=1) : 1266 print(g_help_string) 1267 if(help_level): 1268 print("A full list of options for %s:\n" % ps.label) 1269 for opt in self.valid_opts.olist: 1270 print(" %-20s" % (opt.name )) 1271 if (opt.helpstr != ''): 1272 print(" %-20s %s" % \ 1273 (" use:", opt.helpstr.replace("\n","\n %-20s "%' '))) 1274 if (opt.acceptlist): 1275 print(" %-20s %s" % \ 1276 (" allowed:" , str.join(', ', opt.acceptlist))) 1277 if (opt.deflist): 1278 print(" %-20s %s" % \ 1279 (" default:",str.join(' ', opt.deflist))) 1280 return 1 1281 1282 # remove all the temporary files for epi and anat base names 1283 def cleanup(self, rmold=0): 1284 opt = self.user_opts.find_opt('-epi') 1285 if opt == None : 1286 opt = self.user_opts.find_opt('-dset2') 1287 e = afni_name(opt.parlist[0]) 1288 e.to_afni(new_view=dset_view(e.ppve())) 1289 opt = self.user_opts.find_opt('-anat') 1290 if opt == None : 1291 opt = self.user_opts.find_opt('-dset1') 1292 a = afni_name(opt.parlist[0]) 1293 a.to_afni(new_view=dset_view(a.ppve())) 1294 1295 self.fresh_start( \ 1296 ("%s" % (e.out_prefix())), \ 1297 ("%s" % (a.out_prefix())), rmold = rmold, \ 1298 epipath = ps.output_dir, anatpath = ps.output_dir) 1299 return 1 1300 1301 def version(self): 1302 self.info_msg("align_epi_anat version: %s" % self.align_version) 1303 1304 # copy dataset 1 to dataset 2 1305 # show message and check if dset1 is the same as dset2 1306 # return non-zero error if can not copy 1307 def copy_dset(self, dset1, dset2, message, exec_mode): 1308 self.info_msg(message) 1309 if(dset1.input()==dset2.input()): 1310 print("# copy is not necessary") 1311 return 0 1312# if((os.path.islink(dset1.p())) or (os.path.islink(dset2.p()))): 1313 if(dset1.real_input() == dset2.real_input()): 1314 print("# copy is not necessary") 1315 return 0 1316 ds1 = dset1.real_input() 1317 ds2 = dset2.real_input() 1318 ds1s = ds1.replace('/./','/') 1319 ds2s = ds2.replace('/./','/') 1320 if(ds1s == ds2s): 1321 print("# copy is not necessary - both paths are same") 1322 return 0 1323 print("copying from dataset %s to %s" % (dset1.input(), dset2.input())) 1324 dset2.delete(exec_mode) 1325 com = shell_com( \ 1326 "3dcopy %s %s" % (dset1.input(), dset2.out_prefix()), exec_mode) 1327 com.run() 1328 if ((not dset2.exist())and (exec_mode!='dry_run')): 1329 print("** ERROR: Could not rename %s\n" % dset1.input()) 1330 return 1 1331 return 0 1332 1333 1334 1335## BEGIN script specific functions 1336 def process_input(self): 1337 #Do the default test on all options entered. 1338 #NOTE that default options that take no parameters will not go 1339 #through test, but that is no big deal 1340 for opt in self.user_opts.olist: 1341 if (opt.test() == None): ps.ciao(1) 1342 1343 # skull stripping is on by default for anat/dset1 1344 opt = self.user_opts.find_opt('-anat_has_skull') 1345 if opt != None and opt.parlist[0] == 'no': 1346 ps.skullstrip = 0 1347 else: 1348 ps.skullstrip = 1 1349 ps.skullstrip_method = '3dSkullStrip' 1350 opt = self.user_opts.find_opt('-dset1_strip') 1351 if(opt!=None): 1352 ps.skullstrip_method = opt.parlist[0] 1353 if(ps.skullstrip_method=='None'): 1354 ps.skullstrip = 0 1355 1356 # skull stripping is on by default for epi / dset2 1357 opt = self.user_opts.find_opt('-epi_strip') 1358 if opt != None : 1359 ps.epi_strip_method = opt.parlist[0] 1360 else: 1361 ps.epi_strip_method = '3dSkullStrip' 1362 opt = self.user_opts.find_opt('-dset2_strip') 1363 if(opt!=None): 1364 ps.epi_strip_method = opt.parlist[0] 1365 1366 #Allineate extras 1367 opt = self.user_opts.find_opt('-Allineate_opts') 1368 if opt != None: 1369 ps.AlOpt = str.join(' ', opt.parlist) 1370 else: 1371 ps.AlOpt = '' 1372 1373 # rigid body alignment 1374 opt = self.user_opts.find_opt('-rigid_body') 1375 if opt != None: 1376 ps.AlOpt = ps.AlOpt.replace('-warp aff', '-warp shift_rotate') 1377 1378 # rigid body alignment 1379 opt = self.user_opts.find_opt('-rigid_equiv') 1380 if opt != None: 1381 ps.rigid_equiv = 1 1382 1383 opt = self.user_opts.find_opt('-feature_size') 1384 if opt != None: 1385 featuresize = float(opt.parlist[0]) 1386 blursize = 2.0 * featuresize 1387 aljoin = '-twoblur %f -blok "RHDD(%f)"' % (blursize, featuresize) 1388 ps.AlOpt = "%s %s" % (ps.AlOpt, aljoin) 1389 else: 1390 featuresize = 0.0 1391 1392 opt = self.user_opts.find_opt('-rat_align') 1393 if opt != None: 1394 if featuresize == 0.0 : 1395 featuresize = 0.5 1396 blursize = 1.0 1397 aljoin = '-twoblur %f -blok "RHDD(%f)"' % (blursize, featuresize) 1398 ps.AlOpt = "%s %s" % (ps.AlOpt, aljoin) 1399 ps.skullstrip_opt = "-rat" 1400 1401 #big_move? 1402 opt1 = self.user_opts.find_opt('-big_move') 1403 #giant_move? 1404 opt2 = self.user_opts.find_opt('-giant_move') 1405 #ginormous move? 1406 opt3 = self.user_opts.find_opt('-ginormous_move') 1407 if(not (opt1 or opt2 or opt3)): 1408 ps.AlOpt = "%s -onepass " % ps.AlOpt 1409 else: 1410 # resampling has the potential to cut off data 1411 # the cmass flag is used by 3dAllineate to align the center of mass 1412 # first and then resamples 1413 self.resample_flag = 0 1414 1415 if(opt1): 1416 ps.AlOpt = "%s -twopass " % ps.AlOpt 1417 1418 if(opt2 or opt3): 1419 if featuresize == 0.0 : 1420 fsize = 1 1421 else: 1422 fsize = featuresize 1423 1424 ps.AlOpt = \ 1425 "%s -twobest 11 -twopass -VERB -maxrot 45 -maxshf 40 " \ 1426 "-fineblur %s -source_automask+2" % (ps.AlOpt, fsize) 1427 ps.cmass = "cmass" 1428 ps.giant_move = 1 1429 else : 1430 ps.giant_move = 0 1431 if(opt3): # ginormous option 1432 ps.align_centers = 1 1433 1434 #supersize - allow for large scale changes 1435 opt = self.user_opts.find_opt('-supersize') 1436 if(opt): 1437 ps.AlOpt = \ 1438 "%s -maxscl 1.5" % ps.AlOpt 1439 1440 #giant_move? 1441 opt = self.user_opts.find_opt('-flip_giant') 1442 if(opt): ps.flip_giant = 1; 1443 1444 #add edges 1445 opt = self.user_opts.find_opt('-AddEdge') 1446 if opt == None: 1447 ps.AddEdge = 0 1448 else: 1449 ps.AddEdge = 1 1450 1451 1452 #get anat and epi 1453 1454 opt = self.user_opts.find_opt('-epi') 1455 if opt == None: 1456 opt = self.user_opts.find_opt('-dset2') 1457 if opt == None: 1458 print("** ERROR: Must use -epi or -dset2 option\n") 1459 return 0 1460 if self.cost == '': 1461 self.cost = 'lpa' # make default cost lpa for dset1/2 terminology 1462 self.dset_prep_off() 1463 self.epi_base = "0" # align to 0th sub-brick 1464 self.dset1_generic_name = 'dset1' 1465 self.dset2_generic_name = 'dset2' 1466 1467 e = afni_name(opt.parlist[0]) 1468 ps.epi = e 1469 opt = self.user_opts.find_opt('-anat') 1470 if opt == None: 1471 opt = self.user_opts.find_opt('-dset1') 1472 if opt == None: 1473 print("** ERROR: Must use -anat or -dset1 options\n") 1474 ps.ciao(1) 1475 if self.cost == '': 1476 self.cost = 'lpa' # make default cost lpa for dset1/2 terminology 1477 self.dset_prep_off() # assume no motion correction, slice timing correction,... 1478 self.epi_base = "0" # align to 0th sub-brick 1479 self.dset1_generic_name = 'dset1' 1480 self.dset2_generic_name = 'dset2' 1481 1482 a = ps.anat0 = afni_name(opt.parlist[0]) 1483 1484 # if cost has not been set by using dset1/2 or cost options 1485 if self.cost == '' : 1486 self.cost = 'lpc' # set default cost to lpc 1487 1488# if ps.user_opts.find_opt('-fresh'): 1489# ps.fresh_start(e.out_prefix(), a.out_prefix()) 1490 1491 #epi input 1492 if not e.exist(): 1493 print( "** ERROR: Could not find epi / dset2 dataset\n %s " 1494 % e.input()) 1495 ps.ciao(1) 1496 1497 #anat input 1498 if not a.exist(): 1499 print("** ERROR: Could not find anat / dset1 dataset\n %s " 1500 % a.input()) 1501 ps.ciao(1) 1502 1503 #get 3dTshift options 1504 opt = self.user_opts.find_opt('-tshift_opts') 1505 if opt != None: 1506 ps.tshift_opt = str.join(' ', opt.parlist) 1507 else: 1508 ps.tshift_opt = '-cubic' 1509 1510 #get 3dvolreg options 1511 opt = self.user_opts.find_opt('-volreg_opts') 1512 if opt != None: 1513 ps.reg_opt = str.join(' ', opt.parlist) 1514 else: 1515 ps.reg_opt = '' 1516 1517 #get 3dSkullstrip options 1518 opt = self.user_opts.find_opt('-skullstrip_opts') 1519 if opt != None: 1520 ps.skullstrip_opt = str.join(' ',opt.parlist) 1521 1522 opt = self.user_opts.find_opt('-dset1strip_opts') 1523 if opt != None: 1524 ps.skullstrip_opt = str.join(' ',opt.parlist) 1525 1526 opt = self.user_opts.find_opt('-epistrip_opts') 1527 if opt != None: 1528 ps.epistrip_opt = str.join(' ',opt.parlist) 1529 1530 opt = self.user_opts.find_opt('-dset2strip_opts') 1531 if opt != None: 1532 ps.epistrip_opt = str.join(' ',opt.parlist) 1533 1534 #get epi base type for alignment (specific sub-brick/median/mean) 1535 opt = self.user_opts.find_opt('-epi_base') 1536 if opt == None: 1537 opt = self.user_opts.find_opt('-dset2_base') 1538 if opt == None: 1539 if(self.epi_base==None): 1540 ps.error_msg("Must use -epi_base or -dset2_base options") 1541 ps.ciao(1) 1542 else: 1543 ps.epi_base = opt.parlist[0] 1544 else: 1545 ps.epi_base = opt.parlist[0] 1546 1547 #get volreg_base (matches epi_base by default) 1548 opt = self.user_opts.find_opt('-volreg_base') 1549 if opt != None: 1550 self.volreg_flag = 1 1551 ps.volreg_base = opt.parlist[0] 1552 else: 1553 ps.volreg_base = ps.epi_base 1554 1555# may not need this and only the epi_base parameter instead 1556 #get 3dTstat options 1557 opt = self.user_opts.find_opt('-tstat_opts') 1558 if opt != None: 1559 ps.tstat_opt = str.join(' ', opt.parlist) 1560 else: 1561 ps.tstat_opt = '' 1562 1563 #check for various center of mass options 1564 optc = self.user_opts.find_opt('-cmass') 1565 if optc == None : 1566 #if no cmass option entered, partial coverage? 1567 cmass_opts = 0 1568 opt = self.user_opts.find_opt('-partial_coverage') 1569 if opt != None: 1570 ps.cmass = 'cmass+a' 1571 cmass_opts += 1 1572 opt = self.user_opts.find_opt('-partial_axial') 1573 if opt != None: 1574 ps.cmass = 'cmass+xy' 1575 cmass_opts += 1 1576 opt = self.user_opts.find_opt('-partial_coronal') 1577 if opt != None: 1578 ps.cmass = 'cmass+xz' 1579 cmass_opts += 1 1580 opt = self.user_opts.find_opt('-partial_sagittal') 1581 if opt != None: 1582 ps.cmass = 'cmass+yz' 1583 cmass_opts += 1 1584 if cmass_opts > 1: 1585 self.error_msg("Can only use a single partial_xxx coverage option") 1586 else: 1587 ps.cmass = optc.parlist[0] 1588 1589 #get talairached anatomical dataset 1590 opt = self.user_opts.find_opt('-tlrc_apar') 1591 if opt != None: 1592 anat_tlrc = afni_name(opt.parlist[0]) 1593 ps.tlrc_apar = anat_tlrc 1594 if not anat_tlrc.exist(): 1595 self.error_msg("Could not find anat talairach template dataset\n" \ 1596 " %s " % anat_tlrc.out_prefix()) 1597 else: 1598 self.info_msg("Talairach transformed anatomical: %s" % \ 1599 (anat_tlrc.input())) 1600 if not ps.epi2anat: 1601 self.error_msg("No talairach transformation performed unless aligning epi2anat") 1602 else : 1603 ps.tlrc_apar = "" 1604 opt = self.user_opts.find_opt('-tlrc_epar') 1605 if opt != None: 1606 at = afni_name(opt.parlist[0]) 1607 ps.tlrc_epar = at 1608 if not at.exist(): 1609 self.error_msg("Could not find epi talairach template dataset\n %s " 1610 % at.input()) 1611 else : 1612 ps.tlrc_epar = "" 1613 1614 1615 #get pre-transformation matrix 1616 opt = self.user_opts.find_opt('-pre_matrix') 1617 if opt != None: 1618 ps.pre_matrix = opt.parlist[0] 1619 else : 1620 ps.pre_matrix = "" 1621 1622 #align_centers to get pre-transformation matrix 1623 opt = self.user_opts.find_opt('-align_centers') 1624 if opt != None: # shouldn't happen, defaults to no 1625 if((opt.parlist[0]=='yes') or (opt.parlist[0]=='on')): 1626 ps.align_centers = 1 1627 self.info_msg("turning on align_centers") 1628 if ps.pre_matrix != "" : 1629 self.error_msg( 1630 "Can not use both align_centers and pre-matrix transformations") 1631 1632 #get post-transformation matrix 1633 opt = self.user_opts.find_opt('-post_matrix') 1634 if opt != None: 1635 ps.post_matrix = opt.parlist[0] 1636 else : 1637 ps.post_matrix = "" 1638 1639 # check on the children 1640 ps.child_epis = self.user_opts.find_opt('-child_epi') 1641 if (ps.child_epis == None): 1642 ps.child_epis = self.user_opts.find_opt('-child_dset2') 1643 1644 if ps.child_epis != None: 1645 self.info_msg("-child_epi or child_dset2 option given") 1646 if (not ps.epi2anat) : 1647 if (self.save_Al_in or self.save_tsh or self.save_vr or self.save_rep or 1648 self.save_resample or self.save_epi_ns): 1649 self.info_msg( 1650 "Child epis / dset2 will be processed although not aligning\n" 1651 "epi to anat datasets with -epi2anat / -dset2to1 because a\n" 1652 "save preprocessing option such as -save_vr is selected") 1653 for child_epi_name in ps.child_epis.parlist: 1654 child_epi = afni_name(child_epi_name) 1655 # it's 11:00, do you know where your children are? 1656 if not child_epi.exist(): 1657 self.error_msg("Could not find child epi/dset2 %s\n" 1658 % child_epi.input()) 1659 else: 1660 self.info_msg("Found child epi/dset2 %s\n" 1661 % child_epi.input()) 1662 1663 else : 1664 self.error_msg( 1665 "Child epis / dset2 are not processed unless aligning epi\n" 1666 "to anat datasets with -epi2anat / -dset2to1 or a save\n" 1667 "preprocessing option such as -save_vr is also selected") 1668 ps.child_epis = None 1669 1670 ps.child_anats = self.user_opts.find_opt('-child_anat') 1671 if ps.child_anats == None: 1672 ps.child_anats = self.user_opts.find_opt('-child_dset1') 1673 1674 if ps.child_anats != None: 1675 self.info_msg("-child_anat option given") 1676 if (not ps.anat2epi) : 1677 self.error_msg( 1678 "child anat datasets are not processed unless aligning anat\n" 1679 "to epi or dset1to2 datasets") 1680 ps.child_anats = None 1681 else : 1682 for child_anat_name in ps.child_anats.parlist: 1683 child_anat = afni_name(child_anat_name) 1684 # it's 11:00, do you know where your children are? 1685 if not child_anat.exist(): 1686 self.error_msg("Could not find child anat / dset1\n %s " 1687 % child_anat.input()) 1688 else: 1689 self.info_msg( 1690 "Found child anat / dset1 %s\n" % child_anat.input()) 1691 1692 # output resolution options 1693 # EPI output bounding box and voxel size 1694 min_d = self.min_dim_dset(ps.epi) 1695 mast_dxyz = "" # set default output 1696 mast_dset = "SOURCE" 1697 if ps.giant_move == 1 : 1698 mast_dset = "BASE" 1699 mast_dxyz = '-mast_dxyz %f' % min_d 1700 1701 opt = self.user_opts.find_opt('-master_epi') # epi to anat resolution 1702 if opt == None: 1703 opt = self.user_opts.find_opt('-master_dset2') # dset2to1 resolution 1704 1705 if opt != None: 1706 if(opt.parlist[0]!='MIN_DXYZ'): 1707 if(isFloat(opt.parlist[0])): 1708 min_d = float(opt.parlist[0]) 1709 mast_dxyz = '-mast_dxyz %f' % min_d 1710 else: 1711 mast_dset = opt.parlist[0] 1712 if(opt.parlist[0]!='BASE') and (opt.parlist[0]!='SOURCE'): 1713 mast_dxyz = '' 1714 else: 1715 if(opt.parlist[0]=='SOURCE'): 1716 mast_dxyz = '' 1717 else: # for BASE, resample by default to min.dimension too 1718 mast_dxyz = '-mast_dxyz %f' % min_d 1719 1720 1721 opt = self.user_opts.find_opt('-master_epi_dxyz') 1722 if opt == None: 1723 opt = self.user_opts.find_opt('-master_dset2_dxyz') 1724 if opt != None: 1725 if(isFloat(opt.parlist[0])): 1726 mast_dxyz = '-mast_dxyz %s' % opt.parlist[0] 1727 else: 1728 self.info_msg("****Can not use dxyz setting. Use numbers only***") 1729 self.info_msg("Using default minimum dimension spacing") 1730 1731 # set up 3dAllineate output grid - master dataset for bounding box and 1732 # grid resolution 1733 if((mast_dset=="SOURCE") and ps.align_centers): 1734 self.master_epi_option = "%s" % mast_dxyz 1735 self.master_epi_3dAl_center = 1 1736 else: 1737 self.master_epi_option = \ 1738 "-master %s %s" % (mast_dset, mast_dxyz) 1739 self.master_epi_3dAl_center = 0 1740 self.master_epi_dset = mast_dset 1741 1742 # output resolution options 1743 # EPI output in TLRC dataset's bounding box 1744 # voxel size is minimum dimension of EPI by default 1745 min_d = self.min_dim_dset(ps.epi) 1746 mast_dxyz = '-mast_dxyz %f' % min_d 1747 mast_dset = "BASE" 1748 opt = self.user_opts.find_opt('-master_tlrc') # epi to tlrc resolution 1749 1750 if opt != None: 1751 if(opt.parlist[0]!='MIN_DXYZ'): 1752 if(isFloat(opt.parlist[0])): 1753 min_d = float(opt.parlist[0]) 1754 mast_dxyz = '-mast_dxyz %f' % min_d 1755 else: 1756 mast_dset = opt.parlist[0] 1757 if(opt.parlist[0]!='BASE') and (opt.parlist[0]!='SOURCE'): 1758 mast_dxyz = '' 1759 else: 1760 if(opt.parlist[0]=='SOURCE'): 1761 mast_dxyz = '' 1762 else: # for BASE, resample by default to min.dimension too 1763 mast_dxyz = '-mast_dxyz %f' % min_d 1764 1765 opt = self.user_opts.find_opt('-master_tlrc_dxyz') 1766 if opt != None: 1767 if(isFloat(opt.parlist[0])): 1768 mast_dxyz = '-mast_dxyz %s' % opt.parlist[0] 1769 else: 1770 self.info_msg("****Can not use dxyz setting. Use numbers only***") 1771 self.info_msg("Using default minimum dimension spacing") 1772 1773 # set up 3dAllineate output grid - master dataset for bounding box and 1774 # grid resolution 1775 self.master_tlrc_option = \ 1776 "-master %s %s" % (mast_dset, mast_dxyz) 1777 self.master_tlrc_dset = mast_dset 1778 1779 # anat oblique and 3dAllineate output resolution and bounding box 1780 min_d = self.min_dim_dset(ps.anat0) 1781 mast_dxyz = "" # set default output is same as original 1782 mast_dset = "SOURCE" 1783 if ps.giant_move: 1784 mast_dset = "BASE" 1785 mast_dxyz = '-mast_dxyz %f' % min_d 1786 1787 self.master_anat_option = "-newgrid %f" % min_d 1788# self.info_msg("Spacing for anat to EPI deobliquing is %f mm" % min_d) 1789 1790 opt = self.user_opts.find_opt('-master_anat') 1791 if opt == None: 1792 opt = self.user_opts.find_opt('-master_dset1') 1793 1794 if opt != None: 1795 if(opt.parlist[0]!='MIN_DXYZ'): 1796 if(isFloat(opt.parlist[0])): 1797 min_d = float(opt.parlist[0]) 1798 self.master_anat_option = "-newgrid %f" % min_d 1799 mast_dxyz = '-mast_dxyz %f' % min_d 1800 self.info_msg("Spacing for %s to %s obliquing is %f mm" % \ 1801 (self.dset1_generic_name,self.dset2_generic_name, min_d)) 1802 else: # master is not a voxel size, but a dataset 1803 mast_dset = opt.parlist[0] 1804 # master is a specific dataset name 1805 if(opt.parlist[0]!='BASE') and (opt.parlist[0]!='SOURCE'): 1806 self.master_anat_option = "-gridset %s" % opt.parlist[0] 1807 mast_dxyz = '' 1808 # master is 'SOURCE' or 'BASE' 1809 else: 1810 self.info_msg("****master anat option is %s" 1811 % opt.parlist[0]) 1812 self.info_msg( 1813 "****Can not apply BASE or SOURCE as master " 1814 "for deobliquing. Using default min dim spacing"); 1815 if(opt.parlist[0]=='SOURCE'): 1816 mast_dxyz = '' 1817 else: 1818 mast_dxyz = '-mast_dxyz %f' % min_d 1819 1820 opt = self.user_opts.find_opt('-master_anat_dxyz') 1821 if opt == None: 1822 opt = self.user_opts.find_opt('-master_dset1_dxyz') 1823 if opt != None: 1824 if(isFloat(opt.parlist[0])): 1825 mast_dxyz = '-mast_dxyz %s' % opt.parlist[0] 1826 else: 1827 self.info_msg("****Can not use dxyz setting. Use numbers only***") 1828 self.info_msg("Using default min dim spacing") 1829 1830 # set up 3dAllineate output grid - master dataset for bounding box and 1831 # grid resolution 1832 if((mast_dset=="SOURCE") and ps.align_centers): 1833 self.master_anat_3dAl_option = "%s" % mast_dxyz 1834 self.master_anat_3dAl_center = 1 1835 else: 1836 self.master_anat_3dAl_option = \ 1837 "-master %s %s" % (mast_dset, mast_dxyz) 1838 self.master_anat_3dAl_center = 0 1839 self.master_anat_dset = mast_dset 1840 1841 #get deobliquing options 1842 opt = self.user_opts.find_opt('-deoblique_opts') 1843 if opt != None: 1844 ps.deoblique_opt = str.join(' ',opt.parlist) 1845 else: 1846 ps.deoblique_opt = '' 1847 1848 # user says it's okay to overwrite existing files 1849 opt = self.user_opts.find_opt('-overwrite') 1850 if opt != None: 1851 ps.rewrite = 1 1852 1853 opt = self.user_opts.find_opt('-edge') # use internal edges to drive alignment 1854 if (opt != None): 1855 self.edge = 1 1856 opt = self.user_opts.find_opt('-edge_erodelevel') 1857 if opt == None: 1858 ps.erodelevel = 5 1859 else: 1860 ps.erodelevel = opt.parlist[0] 1861 1862 # the cost function is not as important here because the preprocessing 1863 # steps accomodate for dataset differences - least squares, 1864 # mutual information or local pearson correlation are all good choices 1865 opt = self.user_opts.find_opt('-cost') 1866 if opt == None: self.cost = 'lpa' # local Pearson absolute correlation 1867 if(featuresize==0): # set feature size if not already set 1868 featuresize = 0.5 1869 blursize = 2.0 * featuresize 1870 aljoin = '-twoblur %f -blok "RHDD(%f)"' % (blursize, featuresize) 1871 ps.AlOpt = "%s %s" % (ps.AlOpt, aljoin) 1872 1873 opt = self.user_opts.find_opt('-check_cost') 1874 if opt != None: 1875 1876 self.checkcost = str.join(' ', opt.parlist) 1877 else: 1878 self.checkcost = "" 1879 1880 opt = self.user_opts.find_opt('-multi_cost') 1881 if opt != None: 1882 self.multicost = opt.parlist 1883 else: 1884 self.multicost = [] 1885 self.multicost.append(self.cost) # only the original cost 1886 for mcost in self.multicost: 1887 self.info_msg("Multi-cost is %s" % mcost) 1888 1889 1890 opt = self.user_opts.find_opt('-output_dir') # set alternative output directory 1891 if opt != None: 1892 self.output_dir = opt.parlist[0] 1893 # end with a slash 1894 self.output_dir = "%s/" % os.path.realpath(self.output_dir) 1895 print("# User has selected a new output directory %s" % self.output_dir) 1896 com = shell_com(("mkdir %s" % self.output_dir), self.oexec) 1897 com.run() 1898 print("cd %s" % self.output_dir) 1899 if(not self.dry_run()): 1900 os.chdir(self.output_dir) 1901 else : 1902 self.output_dir = "./" # just the current directory 1903 self.anat_dir = self.output_dir 1904 self.epi_dir = self.output_dir 1905 1906 # all inputs look okay - this goes after all inputs. ########## 1907 return 1 1908 1909 # turn off all preprocessing steps 1910 def prep_off(self) : 1911 self.deoblique_flag = 0 1912 self.tshift_flag = 0 1913 self.volreg_flag = 0 1914 self.resample_flag = 0 1915 self.info_msg("turning off deobliquing tshift, volume registration, resampling") 1916 return 1917 1918 # turn off most preprocessing steps - keep deobliquing though 1919 def dset_prep_off(self) : 1920 self.tshift_flag = 0 1921 self.volreg_flag = 0 1922 self.resample_flag = 0 1923 self.info_msg("turning off tshift, volume registration, resampling") 1924 return 1925 1926 1927 # find smallest dimension of dataset in x,y,z 1928 def min_dim_dset(self, dset=None) : 1929 com = shell_com( \ 1930 "3dAttribute DELTA %s" % dset.input(), ps.oexec,capture=1) 1931 com.run() 1932 if ps.dry_run(): 1933 return (1.234567) 1934 1935 # new purty python way (donated by rick) 1936 min_dx = min([abs(float(com.val(0,i))) for i in range(3)]) 1937 1938 # old non-pythonesque way 1939 if 0: 1940 min_dx = 1000.0 # some high number to start 1941 i = 0 1942 while i < 3 : # find the smallest of the 3 dimensions 1943 dx = abs(float(com.val(0,i))) 1944 if (dx<min_dx): 1945 min_dx = dx 1946 i += 1 1947 1948 1949 if(min_dx==0.0): 1950 min_dx = 1.0 1951 return (min_dx) 1952 1953 1954 # determine if dataset has time shifts in slices 1955 def tshiftable_dset( self, dset=None) : 1956 com = shell_com( \ 1957 "3dAttribute TAXIS_OFFSETS %s" % dset.input(), ps.oexec,capture=1) 1958 com.run() 1959 if ps.dry_run(): 1960 return (1) 1961 if(len(com.so)): status = 1 1962 else: status = 0 1963 return (status) 1964 1965 # determine if dataset is oblique 1966 def oblique_dset( self, dset=None) : 1967 com = shell_com( \ 1968 "3dinfo %s | \grep 'Data Axes Tilt:'|\grep 'Oblique'" % dset.input(),\ 1969 ps.oexec,capture=1) 1970 com.run() 1971 if ps.dry_run(): 1972 return (1) # default for dry run is to assume oblique 1973 if(len(com.so)): 1974# print("length of so is %d" % len(com.so)) 1975# oblstr = com.val(len(com.so),4) 1976# if(oblstr=="Oblique"): 1977 self.info_msg( "Dataset %s is ***oblique****" % dset.input()) 1978 return (1) 1979 self.info_msg( "Dataset %s is not oblique" % dset.input()) 1980 return (0) # if here, then not oblique 1981 1982 1983 # parse 1D file of multiple costs for a particular cost 1984 # the data are arranged like this 1985 ## 3dAllineate -allcostX1D results: 1986 ## ___ ls ___ ___ sp ___ ___ mi ___ ___ crM ___ ___ nmi ___ ___ je ___ ___ hel ___ ___ crA ___ ___ crU ___ ___ lss ___ ___ lpc ___ ___ lpa ___ ___ lpc+ ___ ___ ncd ___ 1987 # 0.778644 0.410080 -0.108398 0.778821 0.989730 10.462562 -0.018096 0.882749 0.862151 0.221356 -0.009732 0.990268 0.522726 0.999801 1988 1989 def get_cost(self, costfilename ,costfunction): 1990 # lpc+ZZ (or lpc+zz) evaluated at end of alignment as lpc 1991 if costfunction.lower() == "lpc+zz" : 1992 costfunction = "lpc" 1993 1994 try: 1995 costfile = open(costfilename,'r') 1996 # ignore first line. That should just say "3dAllineate -allcostX1D results" 1997 costfile.readline() 1998 # read the list of costnames in the second line 1999 costnames = costfile.readline() 2000 # remove the underscores and the leading ## 2001 costnamelist = [(x) for x in costnames.split() if x != "___" and x != "#" and x != "##" ] 2002 # read the list of cost values 2003 costlist = costfile.readline() 2004 # convert into list structure 2005 costs = costlist.split() 2006 costfile.close() 2007 # make dictionary of names and costs 2008 costdict = dict(list(zip(costnamelist, costs))) 2009 2010 # be sure cost function is in the dictionary 2011 if not costfunction in costnamelist: 2012 print("** Error processing cost list file %s" % costfilename) 2013 print(" cost '%s' not found in dict\n" \ 2014 " : %s" % (costfunction, costdict)) 2015 2016 # get cost value from dictionary (error handling if it doesn't exist in list) 2017 costvalue = float(costdict[costfunction]) 2018 return (costvalue) 2019 except: 2020 print("ERROR: error reading cost list file") 2021 2022 # find the cost name and then the corresponding value 2023 # for i, name in enumerate(costnamelist): 2024 # if name == costfunction: 2025 # costvalue = costs[i] 2026 # return costvalue 2027 return (1000) 2028 2029 2030 2031# align the anatomical data to the epi data using 3dAllineate 2032# this is the real meat of the program 2033# note for some of the reasoning: 2034# the output dataset is the result of the alignment of the anatomical to the EPI 2035# This is the default and the preferred output for several reasons despite 2036# its not being standard to the usual processing in the past 2037# First, this does not require the EPI data to be resampled other than volume 2038# registration 2039# For medium to large differences of alignment, the EPI data resampled to its 2040# original grid may lose effective resolution. This is caused by the typical 2041# large difference in slice thickness relative to the EPI x,y voxel size within 2042# slice resolution. 2043# Secondly, the anatomical data is usually higher resolution and relatively 2044# isotropic in voxel dimensions 2045# Thirdly, the anatomical dataset is typically used for structural reference 2046# while the EPI voxel values matter for the analysis. Slight blurring of the 2047# anatomical data is relatively unimportant 2048# Fourthly, the anatomical dataset's higher resolution allows for finer 2049# structural alignment (versus downsampling the anatomical to match the EPI) 2050# One could get around these various issues by using the inverse transform as 2051# done in the epi2anat method and then resampling the output grid to a finer 2052# resolution, but this will not usually be necessary. 2053 2054 def align_anat2epi( self, e=None, a=None, m=None, \ 2055 alopt=" -onepass -weight_frac 1.0 -maxrot 6 " \ 2056 "-maxshf 10 -VERB -warp aff ",\ 2057 suf = "_alnd_epi", costfunction = "lpc"): 2058 #m is the weight brick 2059 self.info_msg( "Aligning %s data to %s data" % \ 2060 (ps.dset1_generic_name, ps.dset2_generic_name )) 2061 # for oblique data or pre and post transformed data, save anat to epi transformation 2062 # matrix in separate temporary 1D file 2063 if((ps.obl_a2e_mat!="") or (ps.pre_matrix!="")) or (ps.edge): 2064 o = a.new("%s_temp%s" % (a.out_prefix(), suf)) 2065 self.anat_mat = "%s%s%s_e2a_only_mat.aff12.1D" % \ 2066 (o.p(), ps.anat0.out_prefix(), suf) 2067 else: 2068 # save transformation matrix with original anatomical name,suf,... 2069 self.anat_mat = "%s%s%s_mat.aff12.1D" % (a.p(), ps.anat0.out_prefix(),suf) 2070 if (ps.anat2epi) or (ps.flip): 2071 o = a.new("%s%s" % (ps.anat0.out_prefix(), suf)) # save the permanent data 2072 else: 2073 o = a.new("__tt_%s%s" % (ps.anat0.out_prefix(), suf)) # save temporary copy 2074 2075 ow = a.new("%s%s_wtal" % (a.out_prefix(), suf)) 2076 of = a.new("%s_flip%s" % (ps.anat0.out_prefix(), suf)) 2077 olr = a.new("__tt_%s_lr%s" % (ps.anat0.out_prefix(), suf)) 2078 2079 if(self.master_anat_dset=='BASE'): 2080 o.view = e.view 2081 ow.view = e.view 2082 of.view = e.view 2083 olr.view = e.view 2084 else: 2085 if(self.master_anat_dset=='SOURCE'): 2086 o.view = a.view 2087 ow.view = a.view 2088 of.view = a.view 2089 olr.view = a.view 2090 else: 2091 manat = afni_name(self.master_anat_dset) 2092 o.view = manat.view 2093 ow.view = manat.view 2094 of.view = manat.view 2095 olr.view = manat.view 2096 2097 anatview = o.view 2098 2099 if (not o.exist() or ps.rewrite or ps.dry_run()): 2100 o.delete(ps.oexec) 2101 ow.delete(ps.oexec) 2102 if m: 2103 wtopt = "-wtprefix %s%s -weight %s" % (ow.p(), ow.out_prefix(), m.input()) 2104 else: 2105 wtopt = "-wtprefix %s%s " % (ow.p(),ow.out_prefix()) 2106 if(ps.cmass==""): 2107 cmass = "" 2108 else: 2109 cmass = "-%s" % ps.cmass 2110 if(ps.checkcost==""): 2111 checkstr = "" 2112 else: 2113 checkstr = "-check %s" % ps.checkcost 2114 com = shell_com( \ 2115 "3dAllineate -%s " # costfunction \ 2116 "%s " # weighting \ 2117 "-source %s " \ 2118 "-prefix %s%s -base %s " \ 2119 "%s " # center of mass options (cmass) \ 2120 "-1Dmatrix_save %s " \ 2121 "%s %s %s " # master grid, other 3dAllineate options \ 2122 # (may be user specified) \ 2123 % (costfunction, wtopt, a.input(), o.p(), o.out_prefix(), \ 2124 e.input(), cmass, self.anat_mat, self.master_anat_3dAl_option, \ 2125 alopt, checkstr ), \ 2126 ps.oexec) 2127 # if this fails, notify the user 2128 if com.run(): 2129 print("** 3dAllineate failure") 2130 return None, None 2131 e2a_mat = self.anat_mat 2132 2133 if (ps.flip): 2134 com = shell_com( \ 2135 "3dLRflip -prefix %s%s -overwrite %s" % \ 2136 (olr.p(), olr.out_prefix(),a.input()), ps.oexec) 2137 com.run() 2138 if((ps.flip_giant) and not(ps.giant_move)): 2139 alopt = \ 2140 "%s -twobest 11 -twopass -VERB -maxrot 45 -maxshf 40 " \ 2141 "-source_automask+2" % (ps.AlOpt) 2142 ps.cmass = "cmass" 2143 # flips don't need weight 2144 if m: 2145 wtopt = "-weight %s" % (m.input()) 2146 else: 2147 wtopt = "" 2148 2149 # save transformation matrix with original anatomical name,suf,... 2150 self.flip_mat = "%s%s_flip_%s_mat.aff12.1D" % \ 2151 (a.p(), ps.anat0.out_prefix(),suf) 2152 2153 com = shell_com( \ 2154 "3dAllineate -%s " # costfunction \ 2155 "%s " # weighting \ 2156 "-source %s " \ 2157 "-prefix %s%s -base %s " \ 2158 "%s " # center of mass options (cmass) \ 2159 "-1Dmatrix_save %s " \ 2160 "%s %s %s " # master grid, other 3dAllineate options \ 2161 # (may be user specified) \ 2162 % (costfunction, wtopt, olr.input(), of.p(), of.out_prefix(), \ 2163 e.input(), cmass, self.flip_mat, self.master_anat_3dAl_option, \ 2164 alopt, checkstr ), \ 2165 ps.oexec) 2166 com.run() 2167 2168 com = shell_com( \ 2169 "3dAllineate -allcostX1D IDENTITY __tt_lr_noflipcosts.1D " \ 2170 "%s " # weighting \ 2171 "-source %s " \ 2172 "-base %s " \ 2173 "%s " # center of mass options (cmass) \ 2174 "%s %s %s " # master grid, other 3dAllineate options \ 2175 % (wtopt, o.input(), e.input(), cmass, \ 2176 self.master_anat_3dAl_option, alopt, checkstr ), \ 2177 ps.oexec) 2178 com.run() 2179 com = shell_com( \ 2180 "3dAllineate -allcostX1D IDENTITY __tt_lr_flipcosts.1D " \ 2181 "%s " # weighting \ 2182 "-source %s " \ 2183 "-base %s " \ 2184 "%s " # center of mass options (cmass) \ 2185 "%s %s %s " # master grid, other 3dAllineate options \ 2186 % (wtopt, of.input(), e.input(), cmass, \ 2187 self.master_anat_3dAl_option, alopt, checkstr ), \ 2188 ps.oexec) 2189 com.run() 2190 2191 noflipcost = self.get_cost("__tt_lr_noflipcosts.1D",costfunction) 2192 print("No flip cost is %f for %s cost function" % (noflipcost, costfunction)) 2193 flipcost = self.get_cost("__tt_lr_flipcosts.1D",costfunction) 2194 print("Flip cost is %f for %s cost function" % (flipcost, costfunction)) 2195 outname = 'aea_checkflip_results.txt' 2196 try: 2197 f = open(outname, 'w') 2198 f.write("flip_cost_orig : %f\n" % noflipcost) 2199 f.write("flip_cost_flipped : %f\n" % flipcost) 2200 f.write("flip_cost_func : %s\n" % costfunction) 2201 f.write("flip_base: %s\n" % e.pv()) 2202 f.write("flip_dset_orig : %s\n" % o.pv()) 2203 f.write("flip_dset_flipped : %s\n" % of.pv()) 2204 if(flipcost < noflipcost): 2205 print("WARNING: ************ flipped data aligns better than original data\n" \ 2206 "Check for left - right flipping in the GUI ************************") 2207 f.write("flip_guess : DO_FLIP\n") 2208 else: 2209 f.write("flip_guess : NO_FLIP\n") 2210 print("Data does not need flipping") 2211 f.close() 2212 except: 2213 print("WARNING: Could not write to flip text file") 2214 2215 # if rigid_equiv, then save only rigid equivalent alignment 2216 # extract with cat_matvec 2217 if(ps.rigid_equiv): 2218 self.info_msg( \ 2219 "Reducing alignment transformation to rigid equivalent") 2220 com = shell_com( \ 2221 "cat_matvec -ONELINE %s -P > __tt_temp.1D ; sleep 1; \ 2222 mv __tt_temp.1D %s" % \ 2223 (self.anat_mat,self.anat_mat), ps.oexec) 2224 com.run(); 2225 2226 # if not doing alignment for anat2epi, just return now, 2227 # and use the xform matrix later 2228 if (not(ps.anat2epi)): 2229 return o, ow 2230 2231 if((ps.obl_a2e_mat!="") or ps.edge ) : 2232 # save the permanent data 2233 o = afni_name("%s%s%s%s" % \ 2234 (ps.anat_dir, ps.anat0.out_prefix(), suf,ps.anat0.view)) 2235 if (not o.exist() or ps.rewrite or ps.dry_run()): 2236 o.delete(ps.oexec) 2237 else: 2238 self.exists_msg(o.input()) 2239 2240 # for oblique and edge data, 2241 # need to apply matrix to skullstripped anat 2242 if(ps.obl_a2e_mat!=""): 2243 # overall transformation A to E is (E2A^-1 PreShift/Oblique) 2244 # 1Dmatrix_apply takes E to A as input so inverse 2245 # E to A = Obl^-1 E2A 2246 obl_mat = "%s -I" % ps.obl_a2e_mat 2247 e2a_mat = "%s%s%s_mat.aff12.1D" % (a.p(), ps.anat0.out_prefix(),suf) 2248 # combine transformations 2249 # earliest transformation is last as input to cat_matvec 2250 com = shell_com( \ 2251 "cat_matvec -ONELINE %s %s > %s" % \ 2252 (self.anat_mat, obl_mat, e2a_mat), ps.oexec) 2253 com.run(); 2254 self.info_msg( \ 2255 "Combining %s to %s and oblique transformations" % \ 2256 (ps.dset1_generic_name, ps.dset2_generic_name )) 2257 2258 else: # just apply the matrix to the original data (edges) 2259 self.info_msg( "Applying transformation to skullstripped %s" % \ 2260 ps.anat_ns0.input()) 2261 e2a_mat = "%s%s%s_mat.aff12.1D" % (a.p(), ps.anat0.out_prefix(),suf) 2262 # input to cat_matvec 2263 com = shell_com( "mv %s %s" % (self.anat_mat, e2a_mat), ps.oexec) 2264 com.run(); 2265 2266 self.info_msg( \ 2267 "renaming e2a transformation matrix is standard matrix (no obliquity)" ) 2268 self.anat_mat = e2a_mat 2269 2270 com = shell_com( \ 2271 "3dAllineate -base %s -1Dmatrix_apply %s " \ 2272 "-prefix %s%s -input %s %s %s" % \ 2273 ( e.input(), e2a_mat, o.p(),o.out_prefix(),\ 2274 ps.anat_ns0.input(),\ 2275 self.master_anat_3dAl_option, alopt ), ps.oexec) 2276 2277 com.run() 2278 else: 2279 self.exists_msg(o.input()) 2280 # process the children 2281 if ps.child_anats != None: 2282 if (ps.anat2epi): 2283 for child_anat_name in ps.child_anats.parlist: 2284 child_anat = afni_name(child_anat_name) 2285 2286 # skip the parent if it's included 2287 if(child_anat.input()==ps.anat0.input()) : 2288 child_anat_out = afni_name("%s%s_child%s%s" % \ 2289 (child_anat.p(),child_anat.prefix,suf,child_anat.view)) 2290 else: 2291 child_anat_out=afni_name("%s%s%s%s" % \ 2292 (child_anat.p(),child_anat.prefix,suf,child_anat.view)) 2293 2294 child_anat_out.view = anatview # child_anat.view 2295 self.info_msg("Processing child %s: %s" % \ 2296 (ps.dset1_generic_name, child_anat.ppv())) 2297 if (ps.rewrite) : 2298 overwritestr = "-overwrite" 2299 child_anat_out.delete(ps.oexec) 2300 else : 2301 overwritestr = "" 2302 com = shell_com( \ 2303 "3dAllineate -base %s -1Dmatrix_apply %s " \ 2304 "-prefix %s%s -input %s %s %s %s" % \ 2305 ( e.input(), e2a_mat, \ 2306 child_anat_out.p(), child_anat_out.out_prefix(), \ 2307 child_anat.input(), \ 2308 self.master_anat_3dAl_option, alopt, overwritestr), \ 2309 ps.oexec) 2310 2311 com.run() 2312 2313 return o, ow 2314 2315# Some notes on this alignment matrix mechanics are sorely needed 2316# because it gets somewhat involved and takes a good deal of time to explain 2317# and to understand all the possible combinations of cases. 2318# Let's first define a series of "spaces" or coordinate systems where 2319# datasets are in alignment with that particular space's original dataset. 2320# 2321# ac = AnatCard (the usual original anatomical dataset without any obliquity 2322# transformation applied, in what we call cardinal or cardinalized space ) 2323# ao = AnatReal (the deobliqued version of the above transformed by its 2324# IJK_TO_DICOM_REAL transformation matrix) 2325# ec = EPICard (similar to AnatCard but the EPI cardinalized dataset, our starting dset) 2326# eo = EPIReal (similar to AnatReal but the deobliqued EPI dataset) 2327# 2328# The usual case we calculate is the transformation from A_ac to A_ec for the anatomical dataset 2329# or the reverse E_ec to E_ac for the epi dataset where the "subscript" signifies the space 2330# Each transformation is computed for a specific space: 2331# A2E = A_ec2E_ec_base = the transformation computed by 3dAllineate to transform the anatomical in 2332# the EPI's cardinal space to be in alignment with the EPI base sub-brick 2333# E_ec_base = the EPI in its own cardinalized space. This is the sub-brick used by 3dAllineate 2334# and the one used for 3dAllineate and the one used as a reference for 3dvolreg 2335# 2336# The "oblique" transformation of the anatomical dataset that transforms the original cardinalized 2337# anatomical dataset to the EPI's cardinalized space is done inside 3dWarp as 2338# Obl = A_ac2A_ec = (T_eo T_ec^-1)^-1 (T_ao T_ac^-1) 2339# 2340 2341# The time-series motion correction volume registration can be described as a transformation 2342# of the EPI sub-bricks to the EPI base sub-brick transformation 2343# VR = E_ec_i2E_ec_base 2344 2345# An additional initial shift or pre-transformation on the original anatomical dataset 2346# (such as produced by @Align_Centers) 2347# In this implementation, we'll treat this transformation as a substitute for the oblique 2348# transformation A_ac2A_ec 2349 2350# One may also choose an additional post-transformation matrix like the one produced by 2351# 3dTagalign usually instead of a tlrc transformation. In this case, the tlrc transformation 2352# would be defined by the post transformation matrix instead of the Warp transformation 2353# in the tlrc_apar dataset's header. 2354 2355# Both the pre and post transformation matrices could potentially be applied 2356# many different ways, but until there's a need for something different, let's use it only as 2357# a substitute for the obliquity and tlrc warp transformation matrices,respectively. 2358 2359# There are the complications that the 1Dmatrix_save and 1Dmatrix_apply use the reverse of what 2360# you might expect, but at least they're consistent. The matrices that come from those options 2361# (used by both 3dvolreg and 3dAllineate) are the transformation of base to source 2362# (output to input). So each of these transformations is the inverse of the transformation 2363# performed or how to get back to where you were. In this program, there are two usages 2364# 3dvolreg saves the VR^-1 transformation (Ebase->Eiec) 2365# 3dAllineate saves the E2A transformation (A2E^-1 or E_ec2A_ac) 2366# Applying also uses the same format of base to source, 2367# so needs to be inverted from source to base 2368# A linear algebra reminder here is the inverse of the product of matrices is equal to 2369# the reverse order with each matrix inverted 2370# (A B)^-1 = B^-1 A^-1 2371# Also the order of operations is first transformation is on right, then proceeds to left 2372# before inverting for 1Dmatrix_apply 2373 2374# So the usual case (-anat2epi without any extra shift) is A_ac to A_ec 2375# and may include oblique transformation 2376# Case 1: A_ac to E_ec 2377# A_ac (A_ac2A_ec) A_ec (A_ec2E_ec_base) A_ec_base = A_ec2E_ec A_ac2A_ec 2378# = A2E Obl 2379# matrix_apply (inverse) = Obl^-1 E2A 2380# E2A is the output of 3dAllineate of A2E 2381 2382# Case 2: E_ec to A_ac (-epi2anat) 2383# may also include obliquity and volume registration 2384# E_ec_i _E_ec_i2E_ec_base_ E_ec_base _A_ec2E_ec^-1_ A_ec _A_ac2A_ec^-1_ A_ac 2385# = A_ac2A_ec^-1 A_ec2E_ec^-1 E_ec_i2E_ec_base 2386# = Obl^-1 A2E^-1 VR 2387# matrix_apply (inverse) = VR^-1 A2E Obl = VR^-1 E2A^-1 Obl 2388# E2A is the output of 3dAllineate of A2E, VR^-1 is output of 3dvolreg 2389 2390# 2391# Case 3: E_ec_i to A_at (-epi2anat -tlrc_apar ...) 2392# talairach transformation of follower EPI to talairach anat) 2393# E_ec_i (E_ec_i2E_ec_base) E_ec_base (A_ec2E_ec^-1) A_ec (A_ac2A_ec^-1) A_ac (A_ac2A_at) A_at 2394# = A_ac2A_at A_ac2A_ec^-1 A_ec2E_ec^-1 E_ec_i2E_ec_base 2395# = Tlrc Obl^-1 A2E^-1 VR 2396# matrix_apply (inverse) = VR^-1 A2E Obl Tlrc^-1 = VR^-1 E2A^-1 Obl Tlrc^-1 2397# E2A is the output of 3dAllineate of A2E, VR^-1 is output of 3dvolreg 2398# 2399 2400 2401 # align the epi to the anatomical but do it using the inverse 2402 # transformation of the alignment of anat to epi 2403 def align_epi2anat( self, e=None, a=None, \ 2404 alopt="",\ 2405 suf = "_alnd_anat"): 2406 # tlrc space output null by default 2407 t = [] 2408 self.info_msg(" Applying alignment for %s to %s" % (ps.dset2_generic_name, ps.dset1_generic_name)) 2409 2410 o = e.new("%s%s" % (self.epi_afniformat.out_prefix(), suf)) 2411 o.path = ps.output_dir 2412# o = afni_name("%s%s" % (self.epi.out_prefix(), suf)) # was e.out_prefix() here 2413 if(self.master_epi_dset == 'SOURCE'): 2414 o.view = "%s" % e.view 2415 if(not(o.view)) : 2416 o.view = dset_view(self.master_epi_dset) 2417# self.info_msg("o.view is SOURCE VIEW %s\n" % o.view) 2418 else: 2419 if(self.master_epi_dset == 'BASE'): 2420 o.view = "%s" % a.view 2421 if(not(o.view)) : 2422 o.view = dset_view(a.ppve()) 2423# self.info_msg("o.view is BASE VIEW %s\n" % o.view) 2424 2425 else: 2426 mepi = afni_name(self.master_epi_dset) 2427 o.view = mepi.view 2428 if(not(o.view)) : 2429 o.view = dset_view(self.master_epi_dset) 2430# self.info_msg("o.view is OTHER VIEW %s\n" % o.view) 2431 2432# self.info_msg("o.view is %s\n" % o.view) 2433 if(not(o.view)) : 2434 o.view = '+orig' 2435 eview = "%s" % o.view 2436# o.view = '+orig' 2437 2438 # allow overwrite in AFNI commands 2439 if (ps.rewrite) : 2440 owrite = "-overwrite" 2441 else: 2442 owrite = "" 2443 2444 if (not o.exist() or ps.rewrite or ps.dry_run()): 2445 o.delete(ps.oexec) 2446 o.view = eview 2447 # anat_mat = "%s%s_mat.aff12.1D" % (ps.anat0.out_prefix(),suf) 2448 epi_mat = "%s%s%s_mat.aff12.1D" % (o.p(), self.epi_afniformat.out_prefix(),suf) 2449 self.info_msg("Inverting %s to %s matrix" % \ 2450 (ps.dset1_generic_name, ps.dset2_generic_name )) 2451 2452 if(ps.obl_a2e_mat!="") : 2453 oblique_mat = "%s" % ps.obl_a2e_mat 2454 else : 2455 oblique_mat = "" 2456 2457 com = shell_com( \ 2458 "cat_matvec -ONELINE %s %s -I > %s" % \ 2459 ( oblique_mat, ps.anat_mat, epi_mat), ps.oexec) 2460 com.run(); 2461 2462 # concatenate volume registration from epi data 2463 if(ps.volreg_flag): 2464 self.info_msg("Concatenating volreg and %s " \ 2465 "to %s transformations" % \ 2466 (ps.dset2_generic_name, ps.dset1_generic_name )) 2467 2468 epi_mat = "%s%s%s_reg_mat.aff12.1D" % (o.p(), self.epi_afniformat.out_prefix(), suf) 2469 com = shell_com( \ 2470 "cat_matvec -ONELINE %s %s -I %s > %s" % \ 2471 (oblique_mat, ps.anat_mat, self.reg_mat, epi_mat), ps.oexec) 2472 com.run(); 2473 2474 2475 self.info_msg( "Applying transformation of %s to %s" % \ 2476 (ps.dset2_generic_name, ps.dset1_generic_name )) 2477 com = shell_com( \ 2478 "3dAllineate -base %s -1Dmatrix_apply %s " \ 2479 "-prefix %s%s -input %s %s %s %s" % \ 2480 ( a.input(), epi_mat, o.p(), o.out_prefix(), e.input(),\ 2481 self.master_epi_option, alopt, owrite), ps.oexec) 2482 2483 com.run() 2484 2485 # mark as not oblique if deobliqued 2486 if(oblique_mat!="") : 2487 o.view = eview 2488 com = shell_com ("3drefit -deoblique %s" % (o.input()), ps.oexec) 2489 com.run() 2490 2491 2492# if (not o.exist() and not ps.dry_run()): 2493# self.error_msg("Could not apply transformation to epi data") 2494# return None 2495 2496 # concatenate talairach or post transformation 2497 if((ps.tlrc_apar != "") or (ps.post_matrix !="")): 2498 2499 self.info_msg( "Concatenating talairach/post, volume registration," \ 2500 " %s to %s transformations" % \ 2501 (ps.dset2_generic_name, ps.dset1_generic_name )) 2502 2503 2504 if(ps.post_matrix != ""): 2505 anat_tlrc_mat = ps.post_matrix 2506 epi_mat = "%s%s%s_post_mat.aff12.1D" % (o.p(),self.epi_afniformat.out_prefix(), suf) 2507 2508 if(ps.tlrc_apar != ""): # tlrc parent trumps post_matrix 2509 com = shell_com( \ 2510 "3dAttribute WARP_TYPE %s" % \ 2511 ps.tlrc_apar.input(), ps.oexec, capture=1) 2512 com.run(); 2513 if(com.status != 0) : 2514 self.error_msg("Warp type not defined for this dataset: %s" % \ 2515 ps.tlrc_apar.input()) 2516 return o,t 2517 2518 tlrc_type = int(com.val(0,0)) 2519 if(tlrc_type != 0) : 2520 self.error_msg("Can not compute transformations for manually" 2521 " talairached data") 2522 return o,t 2523 2524 anat_tlrc_mat = "%s::WARP_DATA" % (ps.tlrc_apar.input()) 2525 epi_mat = "%s%s%s_tlrc_mat.aff12.1D" % (o.p(),self.epi_afniformat.out_prefix(), suf) 2526 2527 # note registration matrix, reg_mat, can be blank and ignored 2528 com = shell_com( \ 2529 "cat_matvec -ONELINE %s -I %s %s -I %s > %s" % \ 2530 (anat_tlrc_mat, oblique_mat, ps.anat_mat, self.reg_mat, \ 2531 epi_mat), ps.oexec) 2532 com.run(); 2533 2534 if(ps.tlrc_apar!=""): 2535 tlrc_dset = afni_name("%s_tlrc%s+tlrc" % (self.epi_afniformat.prefix, suf)) 2536 tlrc_dset.path = o.p() 2537 # tlrc_dset.view = ps.tlrc_apar.view '+tlrc' 2538 if(self.master_tlrc_dset=='SOURCE'): 2539 tlrc_dset.view = e.view 2540 else: 2541 if(self.master_tlrc_dset=='BASE'): 2542 tlrc_dset.view = ps.tlrc_apar.view 2543 else: 2544 mtlrc = afni_name(self.master_tlrc_dset) 2545 tlrc_dset.view = mtlrc.view 2546 2547 if tlrc_dset.exist(): 2548 tlrc_dset.delete(ps.oexec) 2549 atlrcpost = tlrc_dset 2550 self.info_msg( \ 2551 "Applying transformation of %s to %s tlrc parent" % \ 2552 (ps.dset2_generic_name, ps.dset1_generic_name )) 2553 2554 com = shell_com( \ 2555 "3dAllineate -base %s -1Dmatrix_apply %s " \ 2556 "-prefix %s%s -input %s -verb %s %s %s" % \ 2557 ( ps.tlrc_apar.input(), epi_mat, atlrcpost.p(), atlrcpost.prefix,e.input(),\ 2558 ps.master_tlrc_option, alopt, owrite), ps.oexec) 2559 2560 else: 2561 tlrc_orig_dset = afni_name("%s%s_post%s" % (self.epi_afniformat.out_prefix(), suf)) 2562 tlrc_orig_dset = o.p() 2563 tlrc_orig_dset.view = '+orig' 2564 base_dset = a 2565 if(self.master_tlrc_dset=='SOURCE'): 2566 tlrc_orig_dset.view = e.view 2567 else: 2568 if(self.master_tlrc_dset=='BASE'): 2569 tlrc_orig_dset.view = a.view 2570 else: 2571 mtlrc = afni_name(self.master_tlrc_dset) 2572 tlrc_orig_dset.view = mtlrc.view 2573 base_dset = mtlrc 2574 2575 if tlrc_orig_dset.exist(): 2576 tlrc_orig_dset.delete(ps.oexec) 2577 atlrcpost = tlrc_orig_dset 2578 self.info_msg("Applying post transformation matrix to %s" % \ 2579 ps.dset2_generic_name) 2580 2581 com = shell_com( \ 2582 "3dAllineate -base %s -1Dmatrix_apply %s " \ 2583 "-prefix %s%s -input %s -verb %s %s %s" % \ 2584 ( base_dset.input(), epi_mat, atlrcpost.p(), atlrcpost.input(), e.input(),\ 2585 ps.master_tlrc_option, alopt, owrite), ps.oexec) 2586 2587 com.run() 2588 2589 # remove obliquity from output 2590 if(atlrcpost.type == 'BRIK'): 2591 # force +tlrc output for master SOURCE option - 3dAllineate saves this as +orig 2592 if((ps.master_tlrc_dset=="SOURCE") and (ps.tlrc_apar!="")): 2593 com = shell_com ("3drefit -deoblique -view tlrc %s%s+orig" % \ 2594 ( atlrcpost.p(), atlrcpost.prefix), ps.oexec) 2595 com.run() 2596 else: 2597 if(oblique_mat!=""): 2598 com = shell_com ("3drefit -deoblique %s%s+tlrc" % \ 2599 (atlrcpost.p(), atlrcpost.prefix), ps.oexec) 2600 com.run() 2601 t = atlrcpost 2602 else: 2603 self.exists_msg(o.input()) 2604 2605# if (not o.exist() and not ps.dry_run()): 2606# self.error_msg("Could not apply tlrc transformation to epi data") 2607# return None 2608 2609 2610 return o,t 2611 2612 # reduce EPI dataset to a single representative sub-brick 2613 def tstat_epi(self, e=None, tstat_opt="", prefix = "temp_ts" ): 2614 o = afni_name(prefix) 2615 o.to_afni(new_view=dset_view(o.ppve())) 2616 if (not o.exist() or ps.rewrite or ps.dry_run()): 2617 o.delete(ps.oexec) 2618 # if more than 1 sub-brick 2619 if (ps.dry_run() or \ 2620 ((not ps.dry_run() and dset_dims(e.input())[3] > 1))): 2621 # could be: if number choose bucket else use that as stat 2622 # choose a statistic as representative 2623 self.info_msg("Creating representative %s sub-brick" % \ 2624 ps.dset2_generic_name) 2625 2626 # if an integer, choose a single sub-brick 2627 if(ps.epi_base.isdigit()): 2628 # if an integer, choose a single sub-brick 2629 com = shell_com( \ 2630 "3dbucket -prefix %s %s'[%s]'" % \ 2631 (o.pp(), e.input(), ps.epi_base) , ps.oexec) 2632 else: 2633 if((ps.epi_base=='median') or (ps.epi_base=='max') or \ 2634 (ps.epi_base=='mean')): 2635 com = shell_com( \ 2636 "3dTstat -%s -prefix %s %s" % \ 2637 (ps.epi_base, o.pp(), e.input()), ps.oexec) 2638 else: 2639 self.info_msg( 2640 "using 0th sub-brick - assuming epi_base is dataset name" ) 2641 com = shell_com( "3dbucket -prefix %s %s'[0]'" % \ 2642 (o.pp(), e.input()), ps.oexec) 2643 else: # choose a single sub-brick (sub-brick 0) 2644 self.info_msg("using 0th sub-brick because only one found") 2645 com = shell_com( "3dbucket -prefix %s %s'[0]'" % \ 2646 (o.pp(), e.input()), ps.oexec) 2647 com.run(); 2648 if (not o.exist() and not ps.dry_run()): 2649 o.show 2650 print("** ERROR: Could not 3dTstat epi") 2651 return None 2652 else: 2653 self.exists_msg(o.input()) 2654 2655 return o 2656 2657 # create edge dataset of internal brain structure edges 2658 def edge_dset(self, e=None, prefix = "temp_edge", binarize=0,erodelevel=2): 2659 o = afni_name(prefix) 2660 o.view = e.view 2661 m = afni_name("%s_edge_mask" % o.prefix) 2662 m.view = e.view 2663 2664 if (not m.exist() or ps.rewrite or ps.dry_run()): 2665 m.delete(ps.oexec) 2666 self.info_msg("Creating edge dataset") 2667 com = shell_com("3dAutomask -overwrite -erode %s -prefix %s %s" \ 2668 % (erodelevel, m.prefix, e.input()), ps.oexec) 2669 com.run() 2670 2671 if(binarize): 2672 lprefix = "%s_edge_cvar" % o.prefix 2673 else: 2674 lprefix = o.prefix 2675 2676 com = shell_com( \ 2677 "3dLocalstat -overwrite -mask %s" \ 2678 " -nbhd 'RECT(-2,-2,-1)'" \ 2679 " -stat cvar -prefix %s %s" % \ 2680 (m.ppv(), lprefix, e.input()), ps.oexec) 2681 com.run(); 2682 2683 if(binarize): 2684 com = shell_com( \ 2685 "3dhistog -omit 0 -max 1 -nbins 1000 %s_edge_cvar%s " \ 2686 " > %s_edge_histo.1D" % (prefix, o.view, prefix), ps.oexec) 2687 com.run() 2688 2689 com = shell_com("3dTstat -argmax -prefix %s_edge_histomax " \ 2690 "%s_edge_histo.1D'[1]'\\' " % \ 2691 (prefix, prefix), ps.oexec) 2692 com.run() 2693 2694 com = shell_com("1dcat %s_edge_histomax.1D" % \ 2695 prefix, ps.oexec, capture=1) 2696 com.run() 2697 2698 if(ps.dry_run()): edgeindex = 123 2699 else: edgeindex = int(com.val(0,0)) 2700 2701 com = shell_com("1dcat %s_edge_histo.1D'[0]{%d}'" % \ 2702 (prefix, edgeindex), ps.oexec, capture=1) 2703 com.run() 2704 2705 if(ps.dry_run()): edgevalue = 0.0123 2706 else: edgevalue = float(com.val(0,0)) 2707 2708 # threshold anatomical and EPI edge data 2709 self.info_msg("Thresholding edges at %f" % edgevalue) 2710 com = shell_com( \ 2711 '3dcalc -a %s_edge_cvar%s -overwrite -expr "step(a-%f)"' \ 2712 ' -prefix %s' % (prefix, o.view, edgevalue, o.out_prefix()), ps.oexec) 2713 com.run() 2714 2715 if (not o.exist() and not ps.dry_run()): 2716 print("** ERROR: Could not create edge dataset %s" % o.ppv()) 2717 return None 2718 else: 2719 self.exists_msg(o.input()) 2720 2721 2722 return o 2723 2724 2725 # deoblique epi dataset 2726 def deoblique_epi(self, e=None, deoblique_opt="", prefix="temp_deob"): 2727 o = e.new(prefix) 2728 if (not o.exist() or ps.rewrite or ps.dry_run()): 2729 o.delete(ps.oexec) 2730 self.info_msg( "Deobliquing") 2731 com = shell_com( \ 2732 "3dWarp -deoblique -prefix %s%s %s %s " % \ 2733 ( o.p(),o.out_prefix(), deoblique_opt, e.input()), ps.oexec) 2734 com.run(); 2735 if (not o.exist() and not ps.dry_run()): 2736 print("** ERROR: Could not deoblique epi data\n") 2737 return None 2738 else: 2739 self.exists_msg(o.input()) 2740 2741 return o 2742 2743 # oblique anat to epi dataset 2744 # even if neither is really oblique, this shouldn't hurt. 2745 # if a pre-transformation matrix is defined, apply it here 2746 def oblique_anat2epi(self,a=None,e=None,oblique_opt="",suffix="_ob"): 2747 o = a.new("%s%s" % (a.out_prefix(), suffix)) 2748 2749 if(self.pre_matrix!=""): 2750 # use pre-shift transformation matrix instead 2751 self.obl_a2e_mat = self.pre_matrix 2752 warp_str = "3dWarp -verb -matvec_in2out %s -prefix %s%s %s %s %s " \ 2753 % (self.pre_matrix, o.p(), o.out_prefix(), \ 2754 self.master_anat_option, oblique_opt, \ 2755 a.input()) 2756 else: 2757 if(ps.align_centers): 2758 # use shift transformation of centers between grids as initial 2759 # transformation. @Align_Centers (3drefit) instead of 3dWarp 2760 copy_cmd = "3dcopy %s %s%s" % (a.input(), o.p(), o.out_prefix()) 2761 warp_str = "%s; @Align_Centers -base %s -dset %s -no_cp" % \ 2762 (copy_cmd, e.input(), o.input()) 2763 else: 2764 # get obliquity matrix from 3dWarp output and oblique anat to epi 2765 # tempmat = "__tt_%s_obla2e_mat.1D" % a.out_prefix() 2766 self.obl_a2e_mat = "%s%s_obla2e_mat.1D" % (a.p(), a.out_prefix()) 2767 self.info_msg( "Matching obliquity of %s to %s" % \ 2768 (ps.dset1_generic_name, ps.dset2_generic_name )) 2769 2770 warp_str = "3dWarp -verb -card2oblique %s -prefix %s%s %s %s %s " \ 2771 " | \grep -A 4 '# mat44 Obliquity Transformation ::'" \ 2772 " > %s" \ 2773 % (e.input(), o.p(), o.out_prefix(), \ 2774 self.master_anat_option, oblique_opt, \ 2775 a.input(), self.obl_a2e_mat) 2776 2777 if (not o.exist() or ps.rewrite or ps.dry_run()): 2778 o.delete(ps.oexec) 2779 com = shell_com( warp_str, ps.oexec) 2780 com.run(); 2781 if (not o.exist() and not ps.dry_run()): 2782 print("** ERROR: Could not oblique/shift anat to epi data\n") 2783 return None 2784 else: 2785 self.exists_msg(o.input()) 2786 2787 if(ps.align_centers): 2788 # align_centers saves the shift transformation, but 2789 # we use the opposite direction for obliquity, so invert 2790 # the shift transformation too and use that like the 2791 # obliquity matrix in later concatenations 2792 self.obl_a2e_mat = "%s%s_shft_I.1D" % (a.p(), a.out_prefix()) 2793 catcom = "cat_matvec %s%s_shft.1D -I > %s" % \ 2794 (o.p(),o.out_prefix(), self.obl_a2e_mat) 2795 com = shell_com( catcom, ps.oexec) 2796 com.run(); 2797 if(self.master_anat_3dAl_center == 1): 2798 self.info_msg("Using align_centered %s as anatomical master" % \ 2799 o.prefix) 2800 self.master_anat_3dAl_option = "-master %s %s" % \ 2801 (o.input(), self.master_anat_3dAl_option) 2802 if(self.master_epi_3dAl_center == 1): 2803 tempshft = e.new("%s_shft%s" % (e.out_prefix(), suffix)) 2804 2805 copy_cmd = "3dcopy %s %s%s" % \ 2806 (e.input(), tempshft.p(), tempshft.out_prefix()) 2807 warp_str = "%s; @Align_Centers -base %s -dset %s -no_cp" % \ 2808 (copy_cmd, e.input(), tempshft.input()) 2809 2810 if (not tempshft.exist() or ps.rewrite or ps.dry_run()): 2811 tempshft.delete(ps.oexec) 2812 com = shell_com( warp_str, ps.oexec) 2813 com.run(); 2814 if (not tempshft.exist() and not ps.dry_run()): 2815 print("** ERROR: Could not oblique/shift anat to epi data\n") 2816 return None 2817 else: 2818 self.exists_msg(tempshft.input()) 2819 self.info_msg("Using align_centered %s as epi master" % \ 2820 tempshft.prefix) 2821 self.master_epi_option = "-master %s %s" % \ 2822 (tempshft.input(), self.master_epi_option) 2823 2824 return o 2825 2826 2827 # do time shifting of EPI dataset 2828 def tshift_epi( self, e=None, tshift_opt="-cubic", prefix="temp_tsh"): 2829 o = afni_name(prefix) 2830 2831 if (not o.exist() or ps.rewrite or ps.dry_run()): 2832 o.delete(ps.oexec) 2833 self.info_msg( "Correcting for slice timing") 2834 com = shell_com( \ 2835 "3dTshift -prefix %s %s %s " % \ 2836 ( o.pp(), tshift_opt, e.input()), ps.oexec) 2837 com.run(); 2838 if (not o.exist() and not ps.dry_run()): 2839 print("** ERROR: Could not do time shifting of epi data\n") 2840 return e 2841 else: 2842 self.exists_msg(o.input()) 2843 2844 return o 2845 2846 2847 # do volume registration of EPI dataset 2848 def register_epi(self, e=None, reg_opt="-quintic", prefix="temp_vr", \ 2849 motion_prefix = "temp_vr", childflag=0): 2850 o = afni_name(prefix) 2851 2852 if (not o.exist() or ps.rewrite or ps.dry_run()): 2853 o.delete(ps.oexec) 2854 # save the volreg output to file names based on original epi name 2855 # (not temporary __tt_ names) 2856 self.mot_1D = "%s_motion.1D" % (motion_prefix) # motion parameters output 2857 self.reg_mat = "%s%s_mat.aff12.1D" % (o.p(), o.out_prefix()) # volreg transformation matrix 2858 self.info_msg( "Volume registration for %s data" % \ 2859 ps.dset2_generic_name) 2860 2861 # user option for which registration program (3dvolreg,3dWarpDrive,...) 2862 opt = self.user_opts.find_opt('-volreg_method') 2863 if opt != None: 2864 vrcom = opt.parlist[0] 2865 else: 2866 vrcom = '3dvolreg' 2867 2868 if (vrcom == '3dWarpDrive'): 2869 vrcom = '3dWarpDrive -affine_general' 2870 if (vrcom == '3dAllineate'): 2871 vrcom = '3dAllineate -cost lpa -automask -source_automask' 2872 # find base for registration 2873 # could be: if number just use that as base 2874 # if((ps.volreg_base=='median') or (ps.volreg_base=='max') 2875 # or (ps.volreg_base=='mean')): 2876 # choose a statistic as representative 2877 # if an integer, choose a single sub-brick 2878 if(childflag) : # or (vrcom != "3dvolreg") : 2879 base = "%s'[%s]'" % (ps.epi.input(), ps.volreg_base) 2880 elif(ps.volreg_base.isdigit()): 2881 base = "%s" % ps.volreg_base 2882 if(vrcom != '3dvolreg'): 2883 base = "%s'[%s]'" % (ps.epi.input(), ps.volreg_base) 2884 2885 # otherwise median, mean or max 2886 else: 2887 # if more than 1 sub-brick, compute stat, otherwise use 0th 2888 if 0: # (dset_dims(e.ppve())[3] < 2 ): 2889 self.info_msg("Not enough sub-bricks to compute %s" % \ 2890 ps.volreg_base) 2891 base = "0" 2892 else: 2893 ots = e.new("%s_ts_tempalpha" % o.out_prefix()) 2894 base = "%s'[0]'" % ots.input() 2895 if (not ots.exist() or ps.rewrite): 2896 ots.delete(ps.oexec) 2897 2898 # compute stats, volreg, then recompute stats 2899 com = shell_com( \ 2900 "3dTstat -%s -prefix %s%s %s" % \ 2901 (ps.volreg_base, ots.p(), ots.out_prefix(), e.input()), ps.oexec) 2902 com.run() 2903 2904 if(not ots.exist() and not ps.dry_run()): 2905 self.error_msg("Could not create intermediate data" \ 2906 "for time series registration") 2907 ps.ciao(1) 2908 2909 ovr_alpha = e.new("%s_vr_tempalpha" % o.out_prefix()) 2910 2911 if((vrcom != '3dvolreg') and (ps.volreg_base.isdigit())): 2912 base = "%s'[%s]'" % (ps.epi.input(), ps.volreg_base) 2913 com = shell_com( \ 2914 "%s -prefix %s%s -base %s %s %s " % \ 2915 ( vrcom, ovr_alpha.p(), ovr_alpha.out_prefix(), base, \ 2916 reg_opt, e.input()), ps.oexec) 2917 com.run() 2918 2919 ots = e.new("%s_vrt" % o.out_prefix()) 2920 base = "%s'[0]'" % ots.input() 2921 if (not ots.exist() or ps.rewrite): 2922 ots.delete(ps.oexec) 2923 2924 com = shell_com( \ 2925 "3dTstat -%s -prefix %s%s %s" % \ 2926 (ps.volreg_base, ots.p(), ots.out_prefix(), ovr_alpha.input()), \ 2927 ps.oexec) 2928 com.run() 2929 2930 if(not ots.exist() and not ps.dry_run()): 2931 self.error_msg("Could not create intermediate data" \ 2932 "for time series registration") 2933 ps.ciao(1) 2934 2935 com = shell_com( \ 2936 "%s -1Dfile %s -1Dmatrix_save %s " \ 2937 "-prefix %s%s -base %s %s %s " % \ 2938 ( vrcom, self.mot_1D, self.reg_mat, o.p(), o.out_prefix(), base, \ 2939 reg_opt, e.input()), ps.oexec) 2940 com.run() 2941 2942 if (not o.exist() and not ps.dry_run()): 2943 self.error_msg( "Could not do volume registration") 2944 return None 2945 else: 2946 self.exists_msg(o.input()) 2947 2948 return o 2949 2950 # resample EPI data to match higher resolution anatomical data 2951 def resample_epi( self, e=None, resample_opt="", prefix="temp_rs", \ 2952 subbrick=""): 2953 o = afni_name(prefix) 2954 if (not o.exist() or ps.rewrite or ps.dry_run()): 2955 o.delete(ps.oexec) 2956 self.info_msg( "resampling %s to match %s data" % \ 2957 (ps.dset2_generic_name, ps.dset1_generic_name )) 2958 2959 if (subbrick == ""): 2960 sb = "" 2961 else: 2962 if(subbrick.isdigit()): 2963 sb = "[%s]" % subbrick 2964 else: 2965 sb = "[0]" 2966 2967 com = shell_com( \ 2968 "3dresample -master %s -prefix %s -inset %s'%s' -rmode Cu" \ 2969 % (ps.anat_ns.ppv(), o.pp(), e.input(),sb), ps.oexec) 2970 com.run() 2971 if (not o.exist() and not ps.dry_run()): 2972 print("** ERROR: Could not resample\n") 2973 return None 2974 else: 2975 self.exists_msg(o.input()) 2976 2977 return o 2978 2979 # remove skull or outside brain area 2980 def skullstrip_data(self, e=None, use_ss='3dSkullStrip', \ 2981 skullstrip_opt="", prefix = "temp_ns"): 2982 self.info_msg( "removing skull or area outside brain") 2983 if (use_ss == '3dSkullStrip'): #skullstrip epi 2984 n = afni_name(prefix) 2985 if (not n.exist() or ps.rewrite or ps.dry_run()): 2986 n.delete(ps.oexec) 2987 com = shell_com( \ 2988 "3dSkullStrip -orig_vol %s -input %s -prefix %s" \ 2989 % (skullstrip_opt, e.input(), n.pp()) , ps.oexec) 2990 com.run() 2991 if (not n.exist() and not ps.dry_run()): 2992 print("** ERROR: Could not strip skull\n") 2993 return None 2994 else: 2995 self.exists_msg(n.input()) 2996 elif use_ss == '3dAutomask': #Automask epi 2997 n = afni_name(prefix) 2998 j = afni_name("%s__tt_am_%s" % (n.p(),n.pve())) 2999 if (not n.exist() or ps.rewrite or ps.dry_run()): 3000 n.delete(ps.oexec) 3001 com = shell_com( \ 3002 "3dAutomask %s -apply_prefix %s %s" \ 3003 % ( skullstrip_opt, n.pp(), e.input()), ps.oexec) 3004 com.run() 3005 if (not n.exist() and not ps.dry_run()): 3006 print("** ERROR: Could not strip skull with automask\n") 3007 return None 3008 j.delete(ps.oexec) 3009 else: 3010 self.exists_msg(n.input()) 3011 else: 3012 n = e; 3013 return n 3014 3015 # create weighting volume for matching (upper 10% by default) 3016 def create_weight(self, e, sq=1.0, box='no', \ 3017 binit = 'no', perci = -1.0, \ 3018 fati = -1, suf="_wt"): 3019 #e is a preprocessed epi, no skull 3020 # box, bin and fat mask are not used 3021 a = ps.anat_ns 3022 3023 o = afni_name("%s%s%s%s" % (ps.output_dir,e.out_prefix(), suf,e.view)) 3024 if perci < 0: 3025 perci = 90.0; 3026 self.info_msg( "Computing weight mask") 3027 com_str = "3dBrickStat -automask -percentile %f 1 %f %s" \ 3028 % (perci, perci, e.input()) 3029 if(not ps.dry_run()): 3030 com = shell_com( com_str, ps.oexec, capture=1) 3031 com.run() 3032 th = float(com.val(0,1)) 3033 self.info_msg( "Applying threshold of %f on %s" % (th, e.input())) 3034 else: 3035 com = shell_com( com_str, "dry_run") 3036 com.run() 3037 th = -999 3038 self.info_msg( "Would be applying threshold for real run here") 3039 3040 if sq == 1.0: 3041 com = shell_com( 3042 "3dcalc -datum float -prefix %s%s "\ 3043 "-a %s -expr 'min(1,(a/%f))'" \ 3044 % (o.p(), o.out_prefix(), e.input(), th), ps.oexec) 3045 else: 3046 com = shell_com( 3047 "3dcalc -datum float -prefix %s%s -a "\ 3048 "%s -expr 'min(1,(a/%f))**%f'" \ 3049 % (o.p(), o.out_prefix(), e.input(), th, sq), ps.oexec) 3050 3051 if (not o.exist() or ps.rewrite or ps.dry_run()): 3052 o.delete(ps.oexec) 3053 com.run() 3054 else: 3055 self.exists_msg(o.input()) 3056 3057 if (not o.exist() and not ps.dry_run()): 3058 print("** ERROR: Could not create weight dset\n") 3059 return None 3060 3061 return o 3062 3063 3064 # Create double edge images for later viewing 3065 def add_edges(self, d1, d2, d3, d1name, d2name, d3name, listlog = ""): 3066 # specified examinelist for output name or use default in @AddEdge 3067 if (listlog == ""): 3068 aelistlog = "" 3069 else : 3070 aelistlog = "-examinelist %s" % listlog 3071 3072 if (d1name==""): 3073 d1AE = afni_name("%s%s" % (d1.prefix,d1.view)) 3074 else: 3075 d1AE = afni_name("%s%s" % (d1name,d1.view)) 3076 3077 if (d2name==""): 3078 d2AE = afni_name("%s%s" % (d2.prefix,d2.view)) 3079 else: 3080 d2AE = afni_name("%s%s" % (d2name,d2.view)) 3081 3082 if (d3name==""): 3083 d3AE = afni_name("%s%s" % (d3.prefix,d3.view)) 3084 else: 3085 d3AE = afni_name("%s%s" % (d3name,d3.view)) 3086 3087 com = shell_com(("mkdir %sAddEdge"% ps.output_dir), ps.oexec) 3088 com.run() 3089 3090 com = shell_com( \ 3091 "3dcopy -overwrite %s %sAddEdge/%s" % \ 3092 (d1.input(), ps.output_dir, d1AE.prefix), ps.oexec) 3093 com.run() 3094 com = shell_com( \ 3095 "3dcopy -overwrite %s %sAddEdge/%s" % \ 3096 (d2.input(), ps.output_dir, d2AE.prefix), ps.oexec) 3097 com.run() 3098 com = shell_com( \ 3099 "3dcopy -overwrite %s %sAddEdge/%s" % \ 3100 (d3.input(), ps.output_dir, d3AE.prefix), ps.oexec) 3101 com.run() 3102 3103 # changed .input to .shortinput 3 Feb 2012 [rickr,dglen] 3104 com = shell_com( 3105 "cd %sAddEdge; @AddEdge -no_deoblique %s %s %s %s; cd - " \ 3106 % (ps.output_dir,aelistlog, d1AE.shortinput(), d2AE.shortinput(),\ 3107 d3AE.shortinput()), ps.oexec) 3108 com.run() 3109 3110 # do the preprocessing of the EPI data 3111 def process_epi(self, use_ss='3dSkullStrip', childflag=0): 3112 basesuff = "" 3113 if self.epi_dir == '': 3114 basepath = self.epi.p() 3115 else: 3116 basepath = self.epi_dir 3117 3118 if(self.epi.type == 'NIFTI'): 3119 #copy original epi to a temporary file 3120 o = afni_name("%s__tt_%s%s" % \ 3121 (basepath, self.epi.out_prefix(),self.epi.view)) 3122 o.to_afni(new_view=dset_view(self.epi.ppve())) 3123 3124 self.info_msg("Copying NIFTI EPI input to AFNI format") 3125 if (not o.exist() or ps.rewrite or ps.dry_run()): 3126 com = shell_com( "3dcopy %s %s%s" % \ 3127 (self.epi.input(), o.p(), o.pv()), ps.oexec) 3128 com.run(); 3129 basename = o.out_prefix() 3130 baseviewext = "%s%s" % (o.view, o.extension) 3131 basepathname = "%s%s" % (basepath, basename) 3132 tempbasepathname = basepathname 3133 #make filename output without .nii.gz 3134 self.epi_afniformat = afni_name("%s%s%s" % \ 3135 (basepath, self.epi.out_prefix(),self.epi.view)) 3136 self.epi_afniformat.to_afni(new_view=dset_view(self.epi.ppve())) 3137 else: 3138 o = self.epi; 3139 basename = o.out_prefix() 3140 baseviewext = "%s%s" % (o.view, o.extension) 3141 basepathname = "%s%s" % (basepath, basename) 3142 tempbasepathname = "%s__tt_%s" % (basepath, basename) 3143 self.epi_afniformat = self.epi 3144 3145# if Childflag: # no temporary files for children with exceptions below 3146# prepre = "" 3147# else: 3148# prepre = "__tt_" 3149# suff = ps.suffix 3150 3151 # time shift epi data, prepend a prefix 3152 if(self.tshift_flag): 3153 if(self.tshiftable_dset(o)) : 3154 basesuff = "%s_tsh" % basesuff 3155 if(ps.save_tsh): 3156 prefix = "%s%s%s" % (basepathname,basesuff,baseviewext) 3157 else: 3158 prefix = "%s%s%s" \ 3159 % (tempbasepathname,basesuff,baseviewext) 3160 o = self.tshift_epi( o, ps.tshift_opt, prefix=prefix) 3161 else: 3162 self.info_msg("Can not do time shifting of slices. " 3163 "Data is already time shifted") 3164 # if timeshifting was done, this will be the final epi dataset 3165 # and then the concatenated volreg, 3dAllineate transformations will be 3166 # applied 3167 tshift_o = o 3168 # do volume registration 3169 if(self.volreg_flag): 3170 if(ps.dry_run() or \ 3171 (not ps.dry_run() and (dset_dims(o.input())[3] > 1))) : 3172 basesuff = "%s_vr" % basesuff 3173 if(ps.save_vr): 3174 prefix = "%s%s%s%s" % (basepath,self.epi_afniformat.out_prefix(),basesuff,baseviewext) 3175# prefix = "%s%s%s" % (basepathname,basesuff,baseviewext) 3176 else: 3177 prefix = "%s%s%s" \ 3178 % (tempbasepathname,basesuff,baseviewext) 3179 # if aligning epi to anat or saving volreg output, save motion parameters 3180 # if(ps.epi2anat): 3181 motion_prefix = "%s%s%s" % (basepath,self.epi_afniformat.out_prefix(),basesuff) 3182# motion_prefix = "%s%s%s" % (basepath,basename,basesuff) 3183 # else: 3184 # motion_prefix = prefix 3185 o = self.register_epi( o, ps.reg_opt, prefix, motion_prefix, 3186 childflag=childflag) 3187 else: 3188 self.info_msg("Skipping time series volume registration. " 3189 "Must have more than a single sub-brick.") 3190 self.reg_mat = "" # children can be skipped too 3191 3192 volreg_o = o 3193 3194 # if just processing child epi datasets, just go home 3195 # and skip reduction, resampling, skullstripping 3196 if(childflag): 3197 return tshift_o, volreg_o, volreg_o 3198 3199 # reduce epi to a single representative sub-brick 3200 basesuff = "%s_ts" % basesuff 3201 if(ps.save_rep): 3202 prefix = "%s%s%s%s" % (basepath,basename,basesuff,baseviewext) 3203 else: 3204 prefix = "%s%s%s" \ 3205 % (tempbasepathname,basesuff,baseviewext) 3206 3207 o = self.tstat_epi(o, ps.tstat_opt, prefix) 3208 3209 # resample epi to match anat 3210 if(self.resample_flag) : 3211 basesuff = "%s_rs" % basesuff 3212 if(ps.save_resample): 3213 prefix = "%s%s%s%s" % (basepath,basename,basesuff,baseviewext) 3214 else: 3215 prefix = "%s%s%s" \ 3216 % (tempbasepathname,basesuff,baseviewext) 3217 3218 e = self.resample_epi( o,"", prefix) 3219 else: 3220 e = o # no need to resample 3221 3222 # remove outside brain or skull 3223 basesuff = "%s_ns" % basesuff 3224 if(self.save_epi_ns) : 3225 prefix = "%s%s%s" % (basepathname,basesuff,baseviewext) 3226 else: 3227 prefix = "%s%s%s" \ 3228 % (tempbasepathname,basesuff,baseviewext) 3229 3230 skullstrip_o = self.skullstrip_data( e, use_ss, ps.epistrip_opt, prefix) 3231 3232 # use edges to align optionally 3233 if(ps.edge) : 3234 basesuff = "%s_edge" % (basesuff) 3235 if(ps.save_rep): 3236 prefix = "%s%s%s" % (basepathname,basesuff,baseviewext) 3237 else: 3238 prefix = "%s%s%s" % (tempbasepathname,basesuff,baseviewext) 3239 3240 skullstrip_o = self.edge_dset(skullstrip_o, prefix, binarize=0, 3241 erodelevel=ps.erodelevel) 3242 3243 return tshift_o, volreg_o, skullstrip_o 3244 3245 # do the preprocessing of the anatomical data 3246 def process_anat(self): 3247 if self.anat_dir == '': 3248 basepath = self.anat0.p() 3249 else: 3250 basepath = self.anat_dir 3251 3252 self.anat0_master = self.anat0 3253 3254 #copy original anat to a temporary file 3255 if(not ps.flip): 3256 self.anat = afni_name("%s__tt_%s%s" % \ 3257 (basepath, self.anat0.out_prefix(),self.anat0.view)) 3258 else: 3259 self.anat = afni_name("%s%s_unflipped%s" % \ 3260 (basepath, self.anat0.out_prefix(),self.anat0.view)) 3261 3262 self.anat.to_afni(new_view=dset_view(self.anat.ppve())) 3263 3264 if (not self.anat.exist() or ps.rewrite or ps.dry_run()): 3265 com = shell_com( "3dcopy %s %s%s" % \ 3266 (self.anat0.input(), self.anat.p(), self.anat.pv()), ps.oexec) 3267 com.run(); 3268 if (not self.anat.exist() and not ps.dry_run()): 3269 print("** ERROR: Could not copy anat (%d)" % self.anat.exist()) 3270 ps.ciao(1) 3271 else: 3272 self.exists_msg(self.anat.input()) 3273 3274 # now that we have the data in AFNI format, use the AFNI format copy names 3275 self.anat0 = afni_name("%s%s%s" % \ 3276 (basepath, self.anat0.out_prefix(),self.anat0.view)) 3277 3278 self.anat0.to_afni(new_view=dset_view(self.anat.ppve())) 3279 3280 a = self.anat; 3281 3282 #do we need to strip ? 3283 if(ps.skullstrip): 3284 #don't use same type of input. 3285 n = afni_name("%s%s_ns%s" % (a.p(), a.out_prefix(),self.anat.view)) 3286 if (not n.exist() or ps.rewrite or ps.dry_run()): 3287 n.delete(ps.oexec) 3288 self.info_msg( "Removing skull from %s data" % \ 3289 ps.dset1_generic_name) 3290 3291 if(ps.skullstrip_method=="3dSkullStrip"): 3292 com = shell_com( \ 3293 "%s -orig_vol %s -input %s -prefix %s%s" \ 3294 % (ps.skullstrip_method, ps.skullstrip_opt, a.input(), \ 3295 n.p(),n.out_prefix()), ps.oexec) 3296 else: 3297 com = shell_com( \ 3298 "%s %s -apply_prefix %s%s %s" \ 3299 % (ps.skullstrip_method, ps.skullstrip_opt, n.p(), \ 3300 n.out_prefix(), a.input()),ps.oexec) 3301 3302 com.run() 3303 if (not n.exist() and not ps.dry_run()): 3304 print("** ERROR: Could not strip skull\n") 3305 return None 3306 else: 3307 self.exists_msg(o.input()) 3308 else: 3309 n = a 3310 3311 ps.anat_ns0 = n # pre-obliquing or edging skullstripped anat 3312 3313 # match obliquity of anat to epi data 3314 if(self.deoblique_flag or ps.align_centers): 3315 # if either anat or epi is oblique or 3316 # there is a pre-transformation matrix, 3317 # move anat to match epi 3318 if((self.oblique_dset(n)) or (self.oblique_dset(ps.epi)) \ 3319 or (self.pre_matrix!="") or ps.align_centers) : 3320 # set default output spacing if not already set with user options 3321 opt = self.user_opts.find_opt('-master_anat') 3322 if opt == None: 3323 min_d = self.min_dim_dset(ps.anat_ns0) 3324 self.master_3dAl_option = "-mast_dxyz %f" % min_d 3325 self.info_msg( 3326 "Spacing for %s to oblique %s alignment is %f" % 3327 (self.dset1_generic_name, self.dset2_generic_name, min_d)) 3328 3329 n = self.oblique_anat2epi(n, ps.epi, ps.deoblique_opt) 3330 3331 # use edges to match optionally 3332 if(self.edge): 3333 prefix = "%s%s_edge%s" % (a.p(),a.out_prefix(), a.view) 3334 n = self.edge_dset(n, prefix,binarize=0, erodelevel=ps.erodelevel) 3335 3336 #do we need to shift? 3337# optc = self.user_opts.find_opt('-align_centers') 3338# if optc != None and optc.parlist[0] == 'yes': 3339 3340# com = shell_com( \ 3341# "@Align_Centers -cm -base %s -dset %s" \ 3342# % (ps.epi.ppve(), n.ppve() ) , ps.oexec) 3343# com.run() 3344# a_shft = n.new("%s_shft" % n.out_prefix(),"+orig") 3345# if (not a_shft.exist() and not ps.dry_run()): 3346# print "** ERROR: Could not shift anat (%d)" % a_shft.exist() 3347# return None 3348# ps.anat_ns = a_shft 3349# else: 3350 ps.anat_ns = n; 3351 return 1 3352 3353 # create final output files (non-temporary names) 3354 def create_output(self, aae, w, eaa, eaat, suf, epi_in=None, anat_in=None): 3355 3356 #Create a properly named version of anatomy aligned to EPI 3357 opt = ps.user_opts.find_opt('-anat') 3358 if opt == None : 3359 opt = self.user_opts.find_opt('-dset1') 3360 ain = afni_name(opt.parlist[0]) 3361 opt = ps.user_opts.find_opt('-epi') 3362 if opt == None: 3363 opt = ps.user_opts.find_opt('-dset2') 3364 3365 ein = afni_name(opt.parlist[0]) 3366 3367 # save skull stripped anat before alignment 3368 # Yikes! Rick noticed the saved skullstrip was the obliqued one 3369 # and not the original. Should be original ps.anat_ns0, not ps.anat_ns 3370 if(ps.skullstrip and ps.save_skullstrip): 3371 ao_ns = afni_name("%s%s_ns%s" % (ps.anat_ns0.p(),ain.prefix,ain.view)) 3372 self.copy_dset( ps.anat_ns0, ao_ns, 3373 "Creating final output: skullstripped %s data" % \ 3374 self.dset1_generic_name, ps.oexec) 3375 3376 # maybe save pre-obliqued skull stripped anat before alignment 3377 # 3 Aug 2011 [rickr] 3378 if(ps.skullstrip and ps.save_origstrip): 3379 ao_ons = afni_name("%s%s%s%s" % (ps.anat_ns0.p(), ain.prefix, ps.save_origstrip,ain.view)) 3380 self.copy_dset( ps.anat_ns0, ao_ons, 3381 "Creating final output: skullstripped original %s data" % \ 3382 self.dset1_generic_name, ps.oexec) 3383 3384 # save anatomy aligned to epi 3385 if (aae): 3386 # save aligned anatomy 3387 o = afni_name("%s%s%s%s" % (aae.p(),ain.prefix, suf,aae.view)) 3388 o.view = aae.view 3389 self.copy_dset( aae, o, 3390 "Creating final output: %s data aligned to %s" % \ 3391 (self.dset1_generic_name, self.dset2_generic_name) , ps.oexec) 3392 self.save_history(o,ps.oexec) 3393 3394 if (ps.save_tsh): 3395 # save the timeshifted EPI data 3396 if (self.epi_ts and self.tshift_flag) : 3397 eo = afni_name("%s%s_tshft%s" % \ 3398 (self.epi_ts.p(),ein.prefix,self.epi_ts.view)) 3399 self.copy_dset( self.epi_ts, eo, 3400 "Creating final output: time-shifted %s" % \ 3401 self.dset2_generic_name, ps.oexec) 3402 3403 # save the volume registered EPI data 3404 if (ps.save_vr): 3405 if (self.epi_vr and self.volreg_flag): 3406 eo = afni_name("%s%s_vr%s" % 3407 (self.epi_vr.p(),ein.prefix,self.epi_vr.view)) 3408 self.copy_dset( self.epi_vr, eo, 3409 "Creating final output: time series volume-registered %s" % \ 3410 self.dset2_generic_name, ps.oexec) 3411 3412 # save Allineate input datasets 3413 if(ps.save_Al_in): 3414 # save weight used in 3dAllineate 3415 if w: 3416 ow = afni_name("%s%s_wt_in_3dAl%s%s" % \ 3417 (w.p(),ein.prefix,suf,w.view)) 3418 self.copy_dset( w, ow, 3419 "Creating final output: weighting data", ps.oexec) 3420 3421 #save a version of the epi as it went into 3dAllineate 3422 if epi_in: 3423 eo = afni_name("%s%s_epi_in_3dAl%s%s" % \ 3424 (epi_in.p(),ein.prefix, suf,epi_in.view)) 3425 self.copy_dset( epi_in, eo, \ 3426 "Creating final output: " 3427 "%s representative data as used by 3dAllineate" % \ 3428 self.dset2_generic_name, ps.oexec) 3429 3430 #save a version of the anat as it went into 3dAllineate 3431 if anat_in: 3432 ao = afni_name("%s%s_anat_in_3dAl%s%s" % \ 3433 (ps.anat_ns0.p(), ain.prefix, suf, anat_in.view)) 3434 self.copy_dset( anat_in, ao, 3435 "Creating final output: %s data as used by 3dAllineate" % \ 3436 self.dset1_generic_name, ps.oexec) 3437 3438 #Now create a version of the epi that is aligned to the anatomy 3439 if (eaa): 3440 #save the epi aligned to anat 3441 o = afni_name("%s%s%s%s" % \ 3442 (eaa.p(), ein.prefix, suf,eaa.view)) 3443 self.copy_dset( eaa, o, 3444 "Creating final output: %s data aligned to %s" % \ 3445 (self.dset2_generic_name, self.dset1_generic_name), ps.oexec) 3446 self.save_history(o,ps.oexec) 3447 #And a version of the epi that is aligned to the anatomy in standard space 3448 if (eaat): 3449 #save the epi aligned to anat in standard space 3450 #put note in header with command line 3451 self.save_history(eaat,ps.oexec) 3452 3453 return 3454 3455 # remove all the temporary files from a previous run 3456 # remove old results too optionally 3457 def fresh_start(self, epref="", apref="", rmold = 0, epipath="", anatpath="" ): 3458 self.info_msg("Removing all the temporary files") 3459 epref = epref.strip() 3460 apref = apref.strip() 3461 if epref == "" and apref == "": 3462 com = shell_com( "\\rm -f %s__tt_*" % ps.output_dir, ps.oexec) 3463 com.run() 3464 else: 3465 if epref != "": 3466 com = shell_com( "\\rm -f %s__tt_%s*" % (epipath, epref), ps.oexec) 3467 com.run() 3468 if(rmold): 3469 com = shell_com( "\\rm -f %s%s*%s*" % (epipath, epref, ps.suffix), \ 3470 ps.oexec) 3471 com.run() 3472 3473 if apref != "": 3474 com = shell_com( "\\rm -f %s__tt_%s*" % (anatpath, apref), ps.oexec) 3475 com.run() 3476 if(rmold): 3477 com = shell_com( "\\rm -f %s%s*%s*" % (anatpath, apref, ps.suffix), \ 3478 ps.oexec) 3479 com.run() 3480 3481 return 3482 3483 3484 # process children EPI as for parent but 3485 # no alignment of anat to EPI and no representative EPI 3486 # Do time shifting, time series volume registration 3487 def process_child_epi(self, childepi) : 3488 # do preprocessing of epi 3489 # do time shifting, volume registration of child epi data 3490 # if requested; otherwise, just keep the input epi 3491 if(childepi.input()==ps.epi.input()) : # skip the parent if it's included 3492 return # in child list 3493 child = copy.copy(ps) 3494 3495 child.epi = childepi 3496 3497 self.info_msg("Parent %s: Child: %s" % 3498 (ps.epi.input(), child.epi.input())) 3499 3500 child.epi_ts, child.epi_vr, child.epi_ns = \ 3501 child.process_epi(childflag=1) 3502 if (not child.epi_ns): 3503 child.ciao(1) 3504 3505 e = child.epi_ns 3506 a = ps.anat0_master # use parent anat 3507 3508 if(ps.prep_only): # if preprocessing only, exit now 3509 return 3510 3511 if (ps.epi2anat) : # does the user want the epi aligned to the anat 3512 # compute transformation just from applying inverse 3513 child.epi_alnd = \ 3514 child.align_epi2anat(child.epi_ts, a, ps.AlOpt, suf=ps.suffix) 3515 if (not child.epi_alnd): 3516 ps.ciao(1) 3517 else: 3518 child.epi_alnd = '' 3519 3520 3521# Main: 3522if __name__ == '__main__': 3523 3524 3525 ps = RegWrap('align_epi_anat.py') 3526 ps.init_opts() 3527 ps.version() 3528 rv = ps.get_user_opts() 3529 if (rv != None): ps.ciao(1) 3530 3531 #process and check input params 3532 if(not (ps.process_input())): 3533 ps.ciao(1) 3534 3535 # get rid of any previous temporary data 3536 ps.cleanup() 3537 3538 3539 #Now process anatomy and epi 3540 if (not ps.process_anat()): 3541 ps.ciao(1) 3542 # do preprocessing of epi 3543 # final epi2anat option may use timeshifted data as basis 3544 ps.epi_ts, ps.epi_vr, ps.epi_ns = \ 3545 ps.process_epi(use_ss=ps.epi_strip_method) 3546 # (ps.user_opts.find_opt('-epi_strip').parlist[0])) 3547 if (not ps.epi_ns): 3548 ps.ciao(1) 3549 3550 e = ps.epi_ns 3551 ps.epi_alnd_tlrc = None 3552 3553 #Create a weight for final pass 3554 if(not(ps.edge)): 3555 ps.epi_wt = \ 3556 ps.create_weight( e, float(ps.sqmask), ps.boxmask, \ 3557 ps.binmask, ps.perc, -1, suf = "_wt") 3558 else: ps.epi_wt = e # use binary edge image for weight too 3559 3560 if(ps.prep_only): # if preprocessing only, exit now 3561 ps.ciao(0) 3562 3563 #Do alignment to that pesky little epi 3564 for mcost in ps.multicost: 3565 if (mcost==ps.cost): 3566 suff = ps.suffix 3567 else: 3568 suff = "%s_%s" % (ps.suffix, mcost) 3569 3570 a = ps.anat_ns 3571 3572 ps.anat_alnd, ps.anat_alndwt = \ 3573 ps.align_anat2epi(e, a, ps.epi_wt, ps.AlOpt, suff, mcost) 3574 if (not ps.anat_alnd): 3575 ps.ciao(1) 3576 if (ps.epi2anat) : # does the user want the epi aligned to the anat 3577 # compute transformation just from applying inverse 3578 a = ps.anat0_master # especially important for oblique datasets 3579 ps.epi_alnd, ps.epi_alnd_tlrc = \ 3580 ps.align_epi2anat(ps.epi_ts, a, ps.AlOpt, suf=suff) 3581 if (not ps.epi_alnd): 3582 ps.ciao(1) 3583 else: 3584 ps.epi_alnd = '' 3585 3586 if (not ps.anat2epi): 3587 ps.anat_alnd = '' 3588 3589 if(ps.AddEdge): 3590 if(ps.output_dir=="./"): 3591 com = shell_com("\\rm -f AddEdge/*", ps.oexec) 3592 else: 3593 com = shell_com("cd %s; \\rm -f AddEdge/*" % ps.output_dir, ps.oexec) 3594 3595 com.run() 3596 # @AddEdge requires single sub-brick, resampled data (grids must match) 3597 # and skullstripped data to avoid extracranial and cranial edges 3598 # so using data as it went into 3dAllineate here 3599 # check if resampling is turned off, and do that here 3600 # for epi2anat, resample epi_alnd to anat 3601 # if skullstripping is off, we'll have to assume it's already done 3602 3603 # Note 3dAllineate does not require resampling. This script 3604 # provides an option to avoid resampling. 3605 # If @AddEdge does require resampling though, so do it now if it hasn't been done yet 3606 if(ps.resample_flag): 3607 ein_rs = e 3608 else: 3609 baseviewext = "%s%s" % (ps.epi.view, ps.epi.extension) 3610 ein_rs = ps.resample_epi(e, "","%s__tt_%s_rs_in%s" % \ 3611 (ps.output_dir, ps.epi.out_prefix(),baseviewext)) 3612 3613 if (ps.epi2anat and ps.anat2epi): 3614 listlog_a2e = "a2e_examine_list.log" 3615 listlog_e2a = "e2a_examine_list.log" 3616 else : 3617 listlog_a2e = "" 3618 listlog_e2a = "" 3619 3620 3621 if (ps.epi2anat): # if aligning epi to anatomical 3622 # for the edge image, we need a dataset at the resolution of the 3623 # anatomical and the specific epi sub-brick we used as the epi_base 3624 # Use 3dAllineate to create the resampled, skullstripped EPI 3625 # from the same representative sub-brick used to align anat to epi 3626 3627 if((ps.output_dir == "./") or (ps.output_dir == " ")): 3628 outdir = ps.epi.p() 3629 else: 3630 outdir = ps.output_dir 3631 3632 ps.epi_addedge = e.new("__tt_%s_addedge" % ps.epi.out_prefix()) 3633 epi_mat = "%s%s%s_mat.aff12.1D" % \ 3634 (outdir, ps.epi.out_prefix(), ps.suffix) 3635 if((ps.volreg_flag) and (ps.epi_base.isdigit())): 3636 epi_mat = "%s%s%s_reg_mat.aff12.1D{%s}" % \ 3637 (outdir, ps.epi.out_prefix(), suff, ps.epi_base) 3638 3639 # epi_mat = "%s{%s}" % (epi_mat, ps.epi_base) 3640 3641 ps.info_msg( "Applying transformation for %s to %s for @AddEdge" % \ 3642 (ps.dset2_generic_name, ps.dset1_generic_name )) 3643 3644 com = shell_com( \ 3645 "3dAllineate -base %s -1Dmatrix_apply '%s' " \ 3646 "-prefix %s%s -input %s -master BASE" % \ 3647 ( ps.anat0.input(), epi_mat, outdir, \ 3648 ps.epi_addedge.out_prefix(), \ 3649 ein_rs.input()), ps.oexec) 3650 com.run() 3651 ps.add_edges(ps.anat_ns, ein_rs, ps.epi_addedge,\ 3652 "%s_ns" % (ps.anat0.out_prefix()), \ 3653 "%s_ns" % (ps.epi.out_prefix()), \ 3654 "%s%s" % (ps.epi.out_prefix(), suff ), \ 3655 listlog = listlog_e2a) 3656 3657 if (ps.anat2epi): 3658 ps.add_edges(ein_rs, ps.anat_ns, ps.anat_alnd, \ 3659 "%s_ns" % ps.epi.out_prefix(), \ 3660 "%s_ns" % ps.anat0.out_prefix(), \ 3661 "%s%s" % (ps.anat0.out_prefix(), suff), \ 3662 listlog = listlog_a2e) 3663 3664 3665 #Create final results 3666 ps.create_output(ps.anat_alnd, ps.anat_alndwt, ps.epi_alnd, ps.epi_alnd_tlrc, \ 3667 "%s" % suff, e, a) 3668 3669 # process the children 3670 if ps.child_epis != None: 3671 for child_epi_name in ps.child_epis.parlist: 3672 child_epi = afni_name(child_epi_name) 3673 ps.process_child_epi(child_epi) 3674 if (ps.rmrm): # cleanup after the children, remove any extra files 3675 ps.fresh_start(child_epi.out_prefix(), apref="", rmold=0,\ 3676 epipath=("%s" % ps.output_dir) ) 3677 3678 #cleanup after the parents too? 3679 if (ps.rmrm): 3680 ps.cleanup() 3681 3682 3683 print("\n# Finished alignment successfully") 3684 if(ps.AddEdge): 3685 print("To view edges produced by @AddEdge, type:") 3686 print("cd AddEdge") 3687 if (0): #ZSS @AddEdge now launches AFNI 3688 print("afni -niml -yesplugouts &") 3689 3690 if (ps.epi2anat and ps.anat2epi): 3691 print("@AddEdge -examinelist %s" % listlog_a2e) 3692 print("@AddEdge -examinelist %s" % listlog_e2a) 3693 else: 3694 print("@AddEdge") 3695 3696 ps.ciao(0) 3697