1#!/usr/local/bin/python3.8 2 3# python3 status: started 4 5 6# -------------------------------------------------------------------------- 7# [PT: May 30, 2021] update to use all *.nii.gz all the time 8# + so that "AFNI_COMPRESSOR = GZIP" in present functioning 9# doesn't cause probs (though the above will soon not 10# affect *.nii?) 11# [PT: June 2, 2020] ... rolled back to use all *.nii, because RCR fixed how 12# AFNI_COMPRESSOR works with NIFTI (-> now ignores them) 13# -------------------------------------------------------------------------- 14 15 16import sys, os 17import copy 18from time import asctime 19 20# AFNI modules 21from afnipy.afni_base import * 22from afnipy.afni_util import * 23from afnipy.option_list import * 24from afnipy.db_mod import * 25from afnipy import ask_me 26 27g_help_string = """ 28 =========================================================================== 29 auto_warp.py - Nonlinear regisration 30 31 Basic Usage: 32 auto_warp.py -base TT_N27+tlrc -input anat.nii \\ 33 -skull_strip_input yes 34 35 --------------------------------------------- 36 REQUIRED OPTIONS: 37 38 -base : name of reference or template volume 39 -input : name of dataset to be registered 40 41 MAJOR OPTIONS: 42 43 -help : this help message 44 45 OTHER OPTIONS: 46 47 -qblur bB bS : specify 3dQwarp blurs for base and source volumes 48 -qworkhard i0 i1: set the two values for 3dQwarp's -workhard option 49 -qw_opts 'OPTS': Pass all of OPTS as extra options directly to 3dQwarp 50""" 51 52## BEGIN common functions across scripts (loosely of course) 53class RegWrap: 54 def __init__(self, label): 55 self.align_version = "0.06" # software version (update for changes) 56 self.label = label 57 self.valid_opts = None 58 self.user_opts = None 59 self.verb = 1 # a little talkative by default 60 self.save_script = '' # save completed script into given file 61 self.rewrite = 0 #Do not recreate existing volumes 62 self.oexec = "" #dry_run is an option 63 self.anat2temp = 1 # align anat to template by default 64 self.rmrm = 1 # remove temporary files 65 self.prep_only = 0 # do preprocessing only 66 self.odir = os.getcwd() 67 self.resample_flag = 1 # do resample 68 self.deoblique_flag = 1 # deoblique datasets first 69 self.deoblique_opt = "" # deobliquing/obliquing options 70 self.at_opt = "" # @auto_tlrc options 71 self.qw_opts = "" # 3dQwarp extra options 72 self.qworkhard = [0, 1] 73 self.qblur = [-3, -3]# blurs for 3dQwarp 74 self.output_dir = '' # user assigned path for anat and EPI 75 76 return 77 78 79 #Define and initialize defaults of user options 80 def init_opts(self): 81 self.valid_opts = OptionList('init_opts') 82 83 self.valid_opts.add_opt('-base', 1, [], \ 84 helpstr="Template volume.") 85 86 self.valid_opts.add_opt('-input', 1, [], \ 87 helpstr="dataset to be aligned to the template") 88 89 self.valid_opts.add_opt('-keep_rm_files', 0, [], \ 90 helpstr="Don't delete any of the temporary files created here") 91 92 self.valid_opts.add_opt('-prep_only', 0, [], \ 93 helpstr="Do preprocessing steps only without alignment") 94 95 self.valid_opts.add_opt('-help', 0, [], \ 96 helpstr="The main help describing this program with options") 97 98 self.valid_opts.add_opt('-limited_help', 0, [], \ 99 helpstr="The main help without all available options") 100 101 self.valid_opts.add_opt('-option_help', 0, [], \ 102 helpstr="Help for all available options") 103 104 self.valid_opts.add_opt('-version', 0, [], \ 105 helpstr="Show version number and exit") 106 107 self.valid_opts.add_opt('-ver', 0, [], \ 108 helpstr="Show version number and exit") 109 110 self.valid_opts.add_opt('-verb', 1, [], \ 111 helpstr="Be verbose in messages and options" ) 112 113 # 26 Nov 2012 [rickr] 114 self.valid_opts.add_opt('-save_script', 1, [], \ 115 helpstr="save executed script in given file" ) 116 117 self.valid_opts.add_opt('-skip_affine', 1, ['no'], ['yes', 'no'], \ 118 helpstr="Skip the affine registration process\n" \ 119 "Equivalent to -affine_input_xmat ID \n" \ 120 "(apply identity transformation)\n") 121 122 self.valid_opts.add_opt('-skull_strip_base', 1, ['no'], ['yes', 'no'],\ 123 helpstr="Do not skullstrip base/template dataset") 124 self.valid_opts.add_opt('-skull_strip_input', 1, ['no'], ['yes', 'no'],\ 125 helpstr="Do not skullstrip input dataset") 126 127 self.valid_opts.add_opt('-ex_mode', 1, ['script'], \ 128 ['quiet', 'echo', 'dry_run', 'script'], \ 129 helpstr="Command execution mode.\n" \ 130 "quiet: execute commands quietly\n" \ 131 "echo: echo commands executed\n" \ 132 "dry_run: only echo commands\n" ) 133 134 self.valid_opts.add_opt('-overwrite', 0, [],\ 135 helpstr="Overwrite existing files") 136 137 self.valid_opts.add_opt('-suffix', 1,['_al']) 138 139 # child anat datasets 140 self.valid_opts.add_opt('-child_anat', -1,[],\ 141 helpstr="Names of child anatomical datasets") 142 143 144 #resolutions for computing transforms 145 146 # 16 Dec 2013 [rickr] 147 self.valid_opts.add_opt('-qblur', 2, [], 148 helpstr="3dQwarp base and source blurs (FWHM)\n") 149 150 # June 2nd 2014 [ZSS] 151 self.valid_opts.add_opt('-qw_opts', -1, [], 152 helpstr="3dQwarp miscellaneous options.\n" \ 153 "Parameters will get passed directly to 3dQwarp.\n") 154 self.valid_opts.add_opt('-qworkhard', 2, [0, 1], 155 helpstr="3dQwarp -workhard values\n") 156 157 self.valid_opts.add_opt('-warp_dxyz', 1,[0.0],\ 158 helpstr="Resolution used for computing warp (cubic only)\n") 159 self.valid_opts.add_opt('-affine_dxyz', 1,[0.0],\ 160 helpstr="Resolution used for computing initial transform (cubic only)\n") 161 self.valid_opts.add_opt('-affine_input_xmat', 1,deflist=['AUTO'],\ 162 helpstr="Affine transform to put input in standard space.\n"\ 163 "Special values are:\n"\ 164 " 'AUTO' to use @auto_tlrc\n"\ 165 " 'ID' to do nothing\n"\ 166 " 'FILE.1D' for a pre-computed matrix FILE.1D will\n"\ 167 " get applied to the input before Qwarping\n") 168 self.valid_opts.add_opt('-smooth_anat', -1,[],\ 169 helpstr="Smooth anatomy before registration\n") 170 self.valid_opts.add_opt('-smooth_base', -1,[],\ 171 helpstr="Smooth template before registration\n") 172 self.valid_opts.add_opt(name='-unifize_input', npar = 1, deflist='yes',acplist=['yes','no'],\ 173 helpstr="To unifize or not unifize the input\n") 174 self.valid_opts.add_opt('-output_dir', 1,['awpy'],\ 175 helpstr="Set directory for output datasets\n") 176 177 self.valid_opts.add_opt('-followers', npar=-1, deflist=[],\ 178 helpstr="Specify follower datasets\n") 179 self.valid_opts.add_opt('-affine_followers_xmat', npar=-1, deflist=[],\ 180 helpstr="Specify follower datasets' affine transforms\n") 181 182 self.valid_opts.add_opt('-skullstrip_opts', -1, [], 183 helpstr="3dSkullstrip miscellaneous options.\n" \ 184 "Parameters will get passed directly to 3dSkullstrip.\n") 185 186 self.valid_opts.add_opt('-at_opts', -1, [], 187 helpstr="@auto_tlrc miscellaneous options.\n" \ 188 "Parameters will get passed directly to @auto_tlrc.\n") 189 190 self.valid_opts.trailers = 0 # do not allow unknown options 191 192 193 def dry_run(self): 194 if self.oexec != "dry_run": 195 return 0 196 else: 197 return 1 198 199 def apply_initial_opts(self, opt_list): 200 opt1 = opt_list.find_opt('-version') # user only wants version 201 opt2 = opt_list.find_opt('-ver') 202 if ((opt1 != None) or (opt2 != None)): 203 # ps.version() 204 ps.ciao(0) # terminate 205 opt = opt_list.find_opt('-verb') # set and use verb 206 if opt != None: self.verb = int(opt.parlist[0]) 207 208 opt = opt_list.find_opt('-save_script') # save executed script 209 if opt != None: self.save_script = opt.parlist[0] 210 211 opt = opt_list.find_opt('-ex_mode') # set execute mode 212 if opt != None: self.oexec = opt.parlist[0] 213 214 opt = opt_list.find_opt('-keep_rm_files') # keep temp files 215 if opt != None: self.rmrm = 0 216 217 opt = opt_list.find_opt('-prep_only') # preprocessing only 218 if opt != None: self.prep_only = 1 219 220 opt = opt_list.find_opt('-help') # does the user want help? 221 if opt != None: 222 ps.self_help(2) # always give full help now by default 223 ps.ciao(0) # terminate 224 225 opt = opt_list.find_opt('-limited_help') # less help? 226 if opt != None: 227 ps.self_help() 228 ps.ciao(0) # terminate 229 230 opt = opt_list.find_opt('-option_help') # help for options only 231 if opt != None: 232 ps.self_help(1) 233 ps.ciao(0) # terminate 234 235 opt = opt_list.find_opt('-suffix') 236 if opt != None: 237 self.suffix = opt.parlist[0] 238 if((opt=="") or (opt==" ")) : 239 self.error_msg("Cannot have blank suffix") 240 ps.ciao(1); 241 242 # 13 Dec, 2013 [rickr] - for Peter Molfese 243 vals, rv = opt_list.get_type_list(float, '-qblur', length=2) 244 if rv: ps.ciao(1) 245 if vals != None: self.qblur = vals 246 247 248 vals, rv = opt_list.get_type_list(float, '-qworkhard', length=2) 249 if rv: ps.ciao(1) 250 if vals != None: 251 if len(vals) == 2: 252 self.qworkhard = vals 253 else: 254 ps.ciao(1) 255 256 opt = opt_list.find_opt('-qw_opts') 257 istr = ' ' 258 if opt != None: 259 self.qw_opts = ' %s' % \ 260 ' '.join(UTIL.quotize_list(opt.parlist, '\\\n%s '%istr, 1)) 261 262 opt = opt_list.find_opt('-warp_dxyz') 263 if opt != None: 264 self.warp_dxyz = float(opt.parlist[0]) 265 if (self.warp_dxyz != 0.0 and (self.warp_dxyz < 0.01 or self.warp_dxyz > 10.0)): 266 self.error_msg("Bad value for -dxyz of %s" % opt.parlist[0]) 267 ps.ciao(1); 268 else: 269 self.error_msg("This should not happen for pre-defined options"); 270 ps.ciao(1) 271 272 opt = opt_list.find_opt('-affine_dxyz') 273 if opt != None: 274 self.affine_dxyz = float(opt.parlist[0]) 275 if (self.affine_dxyz != 0.0 and (self.affine_dxyz < 0.01 or self.affine_dxyz > 10.0)): 276 self.error_msg("Bad value for -dxyz of %s" % opt.parlist[0]) 277 ps.ciao(1); 278 else: 279 self.error_msg("This should not happen for pre-defined options"); 280 ps.ciao(1) 281 282 opt = opt_list.find_opt('-unifize_input') 283 self.unifize_input = opt_is_yes(opt) 284 285 # any affine transform supplied (defaults to auto_tlrc call) 286 opt = opt_list.find_opt('-affine_input_xmat') 287 self.affine_input_xmat = opt.parlist[0] 288 289 # skip_affine is same as input of IDentity matrix 290 opt = opt_list.find_opt('-skip_affine') 291 if opt_is_yes(opt): 292 self.affine_input_xmat = "ID" 293 294 opt = opt_list.find_opt('-followers') 295 if (opt == None): 296 self.followers=[] 297 else: 298 self.followers = [] 299 for fol in opt.parlist: 300 self.followers.append(afni_name(fol)) 301 opt = opt_list.find_opt('-affine_followers_xmat') 302 if (opt == None): 303 self.affine_followers_xmat=[] 304 else: 305 self.affine_followers_xmat = opt.parlist 306 307 if (len(self.followers)): 308 if (len(self.affine_followers_xmat) and \ 309 len(self.affine_followers_xmat)!=len(self.followers)): 310 error_ex("followers and their transforms don't jive") 311 if (len(self.affine_followers_xmat)==0): 312 for kk in self.followers: 313 self.affine_followers_xmat.append("ID") 314 315 #Parse user options 316 def get_user_opts(self): 317 self.user_opts = read_options(sys.argv, self.valid_opts) 318 if self.user_opts == None: return 1 #bad 319 # no options: apply -help 320 if ( len(self.user_opts.olist) == 0 or \ 321 len(sys.argv) <= 1 ) : 322 ps.self_help() 323 ps.ciao(0) # terminate 324 if self.user_opts.trailers: 325 opt = self.user_opts.find_opt('trailers') 326 if not opt: 327 print("** ERROR: seem to have trailers, but cannot find them!") 328 else: 329 print("** ERROR: have invalid trailing args: %s" % opt.parlist) 330 return 1 # failure 331 332 # apply the user options 333 if self.apply_initial_opts(self.user_opts): return 1 334 335 if self.verb > 3: 336 self.show('------ found options ------ ') 337 338 return 339 340 def show(self, mesg=""): 341 print('%s: %s' % (mesg, self.label)) 342 if self.verb > 2: self.valid_opts.show('valid_opts: ') 343 self.user_opts.show('user_opts: ') 344 345 def info_msg(self, mesg=""): 346 if(self.verb >= 1) : 347 print("#++ %s" % mesg) 348 349 def error_msg(self, mesg=""): 350 print("#**ERROR %s" % mesg) 351 352 def error_ex(self, mesg=""): 353 print("#**ERROR %s" % mesg) 354 self.ciao(1) 355 356 def exists_msg(self, dsetname=""): 357 self.info_msg(mesg="** Dataset: %s already exists: Reusing." % dsetname); 358 359 360 def ciao(self, i): 361 if i > 0: 362 print("** ERROR - script failed") 363 elif i==0: 364 print("") 365 366 os.chdir(self.odir) 367 368 if self.save_script: 369 write_afni_com_history(self.save_script) 370 371 sys.exit() 372 373 # save the script command arguments to the dataset history 374 def save_history(self, dset, exec_mode): 375 self.info_msg("Saving history") # sounds dramatic, doesn't it? 376 cmdline = args_as_command(sys.argv, \ 377 '3dNotes -h "', '" %s' % dset.input()) 378 com = shell_com( "%s\n" % cmdline, exec_mode) 379 com.run() 380 381 # show help 382 # if help_level is 1, then show options help only 383 # if help_level is 2, then show main help and options help 384 def self_help(self, help_level=0): 385 if(help_level!=1) : 386 print(g_help_string) 387 if(help_level): 388 print("A full list of options for %s:\n" % ps.label) 389 for opt in self.valid_opts.olist: 390 print(" %-20s" % (opt.name )) 391 if (opt.helpstr != ''): 392 print(" %-20s %s" % \ 393 (" use:", opt.helpstr.replace("\n","\n %-20s "%' '))) 394 if (opt.acceptlist): 395 print(" %-20s %s" % \ 396 (" allowed:" , ', '.join(opt.acceptlist))) 397 if (opt.deflist): 398 if type(opt.deflist[0]) != str: # 31 Mar 2014 [rickr] 399 print(" %-20s %s" % (" default:",opt.deflist)) 400 else: 401 print(" %-20s %s" % \ 402 (" default:",' '.join(opt.deflist))) 403 return 1 404 405 # remove all the temporary files for epi and anat base names 406 def cleanup(self, rmold=0): 407 opt = self.user_opts.find_opt('-output_dir') 408 return 0 409 410 def version(self): 411 self.info_msg("auto_warp.py version: %s" % self.align_version) 412 413 # copy dataset 1 to dataset 2 414 # show message and check if dset1 is the same as dset2 415 # return non-zero error if cannot copy 416 def copy_dset(self, a, prefix): 417 n = afni_name(prefix) 418 if (not n.exist() or ps.rewrite or ps.dry_run()): 419 n.delete(ps.oexec) 420 com = shell_com( \ 421 "3dcopy %s %s" \ 422 % (a.input(), n.input()) , ps.oexec) 423 com.run() 424 if (not n.exist() and not ps.dry_run()): 425 print("** ERROR: Could not copy dset\n") 426 ps.ciao(1) 427 else: 428 self.exists_msg(n.input()) 429 430 return n 431 432 433## BEGIN script specific functions 434 def process_input(self): 435 for opt in self.user_opts.olist: 436 if (opt.test() == None): ps.ciao(1) 437 438 # do not allow AFNI_COMPRESSOR (or avoid NIFTI) 10 Jun 2015 [rickr] */ 439 ename = 'AFNI_COMPRESSOR' 440 if ename in os.environ: 441 print('-- clearing %s ...' % ename) 442 del os.environ[ename] 443 444 opt = self.user_opts.find_opt('-skull_strip_input') 445 if opt == None: 446 self.error_ex("should not be empty") 447 448 if opt.parlist[0] == 'no': 449 ps.skullstrip_input = 0 450 else: 451 ps.skullstrip_input = 1 452 453 454 opt = self.user_opts.find_opt('-skull_strip_base') 455 if opt == None: 456 self.error_ex("should not be empty") 457 458 if opt.parlist[0] == 'no': 459 ps.skullstrip_base = 0 460 else: 461 ps.skullstrip_base = 1 462 463 #get anat and epi 464 465 opt = self.user_opts.find_opt('-input') 466 if opt == None: 467 self.error_ex("No -input"); 468 469 ps.input = afni_name(opt.parlist[0]) 470 if (not ps.input.exist()): 471 self.error_ex("Could not find input") 472 473 opt = self.user_opts.find_opt('-base') 474 if opt == None: 475 self.error_ex("No -template"); 476 477 ps.base = afni_name(opt.parlist[0]) 478 ps.base.locate() 479 if (not ps.base.exist()): 480 self.error_ex("Could not find base") 481 482 #get 3dSkullstrip options 483 opt = self.user_opts.find_opt('-skullstrip_opts') 484 if opt != None: 485 ps.skullstrip_opt = ' '.join(opt.parlist) 486 else: 487 ps.skullstrip_opt = '' 488 489 #get auto_tlrc options 490 opt = self.user_opts.find_opt('-at_opts') 491 if opt != None: 492 ps.at_opt = ' '.join(opt.parlist) 493 else: 494 ps.at_opt = '' 495 496 497 # user says it's okay to overwrite existing files 498 opt = self.user_opts.find_opt('-overwrite') 499 if opt != None: 500 ps.rewrite = 1 501 502 opt = self.user_opts.find_opt('-output_dir') # set alternative output directory 503 if opt != None: 504 self.output_dir = opt.parlist[0] 505 self.output_dir = "%s/" % os.path.realpath(self.output_dir) 506 print("# Output directory %s" % self.output_dir) 507 508 com = shell_com(("mkdir %s" % self.output_dir), self.oexec) 509 com.run() 510 print("cd %s" % self.output_dir) 511 if(not self.dry_run()): 512 os.chdir(self.output_dir) 513 514 # all inputs look okay - this goes after all inputs. ########## 515 return 1 516 517 518 # find smallest dimension of dataset in x,y,z 519 def min_dim_dset(self, dset=None) : 520 com = shell_com( \ 521 "3dAttribute DELTA %s" % dset.input(), ps.oexec,capture=1) 522 com.run() 523 if ps.dry_run(): 524 return (1.234567) 525 526 # new purty python way (donated by rick) 527 min_dx = min([abs(float(com.val(0,i))) for i in range(3)]) 528 529 if(min_dx==0.0): 530 min_dx = 1.0 531 return (min_dx) 532 533 534 # resample EPI data to match higher resolution anatomical data 535 def resample_epi( self, e=None, resample_opt="", prefix="temp_rs", \ 536 subbrick=""): 537 o = afni_name(prefix) 538 if (not o.exist() or ps.rewrite or ps.dry_run()): 539 o.delete(ps.oexec) 540 self.info_msg( "resampling %s to match %s data" % \ 541 (ps.dset2_generic_name, ps.dset1_generic_name )) 542 543 if (subbrick == ""): 544 sb = "" 545 else: 546 if(subbrick.isdigit()): 547 sb = "[%s]" % subbrick 548 else: 549 sb = "[0]" 550 551 com = shell_com( \ 552 "3dresample -master %s -prefix %s -inset %s'%s' -rmode Cu" \ 553 % (ps.anat_ns.ppv(), o.pp(), e.input(),sb), ps.oexec) 554 com.run() 555 if (not o.exist() and not ps.dry_run()): 556 print("** ERROR: Could not resample\n") 557 ps.ciao(1) 558 else: 559 self.exists_msg(o.pve()) 560 561 return o 562 563 # remove skull or outside brain area 564 def skullstrip_data(self, e=None, use_ss='3dSkullStrip', \ 565 skullstrip_opt="", prefix = "temp_ns"): 566 self.info_msg( "removing skull or area outside brain") 567 if (use_ss == '3dSkullStrip'): #skullstrip epi 568 n = afni_name(prefix) 569 if (not n.exist() or ps.rewrite or ps.dry_run()): 570 n.delete(ps.oexec) 571 com = shell_com( \ 572 "3dSkullStrip -orig_vol %s -input %s -prefix %s" \ 573 % (skullstrip_opt, e.input(), n.input()) , ps.oexec) 574 com.run() 575 if (not n.exist() and not ps.dry_run()): 576 print("** ERROR: Could not strip skull\n") 577 ps.ciao(1) 578 else: 579 self.exists_msg(n.input()) 580 elif use_ss == '3dAutomask': #Automask epi 581 n = afni_name(prefix) 582 j = afni_name("%s__tt_am_%s" % (n.p(),n.pve())) 583 if (not n.exist() or ps.rewrite or ps.dry_run()): 584 n.delete(ps.oexec) 585 com = shell_com( \ 586 "3dAutomask -prefix %s %s && 3dcalc -a %s "\ 587 "-b %s -prefix %s -expr 'step(b)*a'" \ 588 % ( j.pp(), e.input(), e.input(), 589 j.input(), n.pp()), ps.oexec) 590 com.run() 591 if (not n.exist() and not ps.dry_run()): 592 print("** ERROR: Could not strip skull with automask\n") 593 ps.ciao(1) 594 j.delete(ps.oexec) 595 else: 596 self.exists_msg(n.input()) 597 else: 598 n = e; 599 return n 600 601 #unifize? 602 def unifize_data(self, a, prefix = "temp_un"): 603 n = afni_name(prefix) 604 if (not n.exist() or ps.rewrite or ps.dry_run()): 605 n.delete(ps.oexec) 606 com = shell_com( \ 607 "3dUnifize -GM -input %s -prefix %s" \ 608 % (a.input(), n.input()) , ps.oexec) 609 com.run() 610 if (not n.exist() and not ps.dry_run()): 611 print("** ERROR: Could not strip skull\n") 612 ps.ciao(1) 613 else: 614 self.exists_msg(n.input()) 615 return n 616 617 # prepare input 618 def prepare_input(self, use_ss='3dSkullStrip'): 619 prefix = "%s/anat.ns.nii" % (ps.output_dir) 620 if (ps.skullstrip_input): 621 skullstrip_o = self.skullstrip_data(self.input, use_ss, ps.skullstrip_opt, prefix) 622 else: 623 skullstrip_o = self.copy_dset (self.input, prefix="%s/anat.nii" % (ps.output_dir)) 624 if (ps.unifize_input): 625 prefix = "%s/anat.un.nii" % (ps.output_dir) 626 skullstrip_o = self.unifize_data(skullstrip_o, prefix) 627 628 return skullstrip_o 629 630 # do the preprocessing of the anatomical data 631 def prepare_base(self, use_ss='3dSkullStrip'): 632 prefix = "%s/base.ns.nii" % (ps.output_dir) 633 if (ps.skullstrip_base): 634 skullstrip_o = self.skullstrip_data( self.base, use_ss, ps.skullstrip_opt, prefix) 635 else: 636 skullstrip_o = self.copy_dset (self.base, prefix="%s/base.nii" % (ps.output_dir)) 637 638 return skullstrip_o 639 640 def resample(self,a,prefix='resampled', dxyz=0.0, m=None): 641 n = afni_name(prefix) 642 if (m == None): 643 m = a 644 if (not n.exist() or ps.rewrite or ps.dry_run()): 645 n.delete(ps.oexec) 646 if (dxyz != 0.0): 647 com = shell_com( \ 648 "3dresample -inset %s -prefix %s -dxyz %f %f %f "\ 649 " -rmode Li -master %s" \ 650 % ( a.input(), n.input(), dxyz, dxyz, dxyz, 651 m.input()), ps.oexec) 652 else: 653 com = shell_com( \ 654 "3dresample -inset %s -prefix %s "\ 655 " -rmode Li -master %s" \ 656 % ( a.input(), n.input(), 657 m.input()), ps.oexec) 658 com.run() 659 if (not n.exist() and not ps.dry_run()): 660 print("** ERROR: Could not strip skull with automask\n") 661 ps.ciao(1) 662 else: 663 self.exists_msg(n.input()) 664 return(n) 665 666 667 def match_resolutions(self, a, b, suf, dxyz=0.0, m=None): 668 if (dxyz != 0.0): 669 print("%f" % (dxyz)) 670 ar = self.resample(a,prefix="anat%s.nii" % suf, dxyz=dxyz, m=m) 671 br = self.resample(b,prefix="base%s.nii" % suf, dxyz=dxyz, m=m) 672 else: 673 ar = a 674 br = b 675 dxa = self.min_dim_dset(a) 676 dxb = self.min_dim_dset(b) 677 print("%f %f" % (dxa, dxb)) 678 if (dxb < dxa): 679 br = self.resample(b,prefix="base%s.nii" % suf, dxyz=dxa, m=m) 680 if (dxa > dxb): 681 ar = self.resample(a,prefix="anat%s.nii" % suf, dxyz=dxb, m=m) 682 683 com = shell_com(\ 684 "3dinfo -same_grid %s %s" % (ar.input(), br.input()), ps.oexec, capture=1) 685 com.run() 686 if (len(com.so) and int(com.so[0]) == 0): 687 ar = self.resample(ar,prefix="anat%sb.nii" % suf ,m=br) 688 689 return ar,br 690 691 692 def align_auto_tlrc( self, a=None, b=None, \ 693 alopt=" ",\ 694 suf = ".aff", at_opt = "", 695 xmat=None): 696 #m is the weight brick 697 self.info_msg( "Aligning %s data to %s data" % \ 698 (b.input(), a.input() )) 699 n = a.new(new_pref="%s%s" % (a.prefix, suf), parse_pref=1) 700 if (not n.exist() or ps.rewrite or ps.dry_run()): 701 n.delete(ps.oexec) 702 if (xmat==None): 703 com = shell_com( \ 704 "@auto_tlrc -base %s " \ 705 " -input %s " \ 706 " -suffix %s " \ 707 " -no_ss -no_pre" \ 708 " -init_xform CENTER %s" \ 709 % ( b.input(), a.input(), suf, at_opt ), \ 710 ps.oexec) 711 else: 712 com = shell_com( \ 713 "3dWarp -matvec_out2in %s -prefix %s %s " \ 714 % ( xmat, n.input(), a.input() ), \ 715 ps.oexec) 716 com.run() 717 if (not n.exist() and not ps.dry_run()): 718 self.error_msg("Failed in affine step"); 719 ps.ciao(1) 720 xmat = "%s.Xat.1D" % n.prefix 721 if (not os.path.isfile(xmat) and not ps.dry_run()): 722 self.error_msg("Failed to find xmat %s" % xmat); 723 ps.ciao(1) 724 else: 725 self.exists_msg(n.input()) 726 if (xmat==None): 727 xmat = "%s.Xat.1D" % n.prefix 728 729 return n,xmat 730 731 def qwarping(self, a=None, b = None, prefix=None): 732 self.info_msg( "Aligning %s data to %s data" % \ 733 (b.input(), a.input() )) 734 if (prefix==None): 735 prefix = "%s.qw" % a.prefix 736 n = a.new(new_pref=prefix, parse_pref=1) 737 if (not n.exist() or ps.rewrite or ps.dry_run()): 738 n.delete(ps.oexec) 739 if self.qworkhard[0] < 0: 740 whopt = '-workhard' 741 else: whopt = "-workhard:%d:%d" % (self.qworkhard[0], self.qworkhard[1]) 742 com = shell_com( \ 743 "3dQwarp "\ 744 " -prefix %s "\ 745 " -blur %s %s %s "\ 746 " %s -base %s -source %s "\ 747 % ( n.input(), self.qblur[0], self.qblur[1], 748 whopt, 749 self.qw_opts, b.input(), a.input()), ps.oexec) 750 com.run() 751 if (not n.exist() and not ps.dry_run()): 752 self.error_msg("Failed in warping step"); 753 ps.ciao(1) 754 else: 755 self.exists_msg(n.input()) 756 757 w = n.new(new_pref="%s_WARP" % n.prefix) 758 759 return (n,w) 760 761 def qwarp_applying(self, a, aff, wrp, prefix=None, dxyz=0.0, master=None): 762 self.info_msg( "Applying warps to %s" % \ 763 ( a.input() )) 764 if (prefix==None): 765 prefix = "./%s.aw.nii" % a.prefix 766 n = afni_name(prefix) 767 n.new_path("./") 768 if (not n.exist() or ps.rewrite or ps.dry_run()): 769 n.delete(ps.oexec) 770 if (aff != "ID"): 771 # waff = "-affter %s" % aff 772 nwarp = '"%s %s"' % (wrp.input(), aff) # 7 Nov, 2014 [rickr] 773 else: 774 # waff = "" 775 nwarp = '%s' % wrp.input() 776 if (dxyz==0.0): 777 dxopt = "" 778 else: 779 dxopt = "-dxyz %f" % dxyz 780 781 # warp datasets may need to grow, so allow for a different master 782 if master == None: mast_str = "NWARP" 783 else: mast_str = "%s" % master.input() 784 785 com = shell_com( \ 786 "3dNwarpApply "\ 787 "-nwarp %s "\ 788 "-master %s "\ 789 " %s "\ 790 "-source %s "\ 791 "-prefix %s "\ 792 % ( nwarp, mast_str, dxopt, a.input(), n.input()), 793 ps.oexec) 794 com.run() 795 if (not n.exist() and not ps.dry_run()): 796 n.show() 797 self.error_msg("Failed in warping step"); 798 ps.ciao(1) 799 else: 800 self.exists_msg(n.input()) 801 return n 802 803 def align_epi_anat(self, e, a, aff): 804 prefix = "%s_al" % e.prefix 805 n = e.new(new_pref=prefix) 806 n.new_path("./") 807 affout = "%s_al_post_mat.aff12.1D" % e.prefix 808 if (not os.path.isfile(affout) or not n.exist() or ps.rewrite or ps.dry_run()): 809 n.delete(ps.oexec) 810 if (aff != "ID"): 811 waff = "-post_matrix %s" % aff 812 else: 813 waff = "" 814 com = shell_com( \ 815 "align_epi_anat.py "\ 816 " -epi %s -epi_base 3 -epi2anat "\ 817 " -anat %s -big_move "\ 818 " -anat_has_skull no -overwrite "\ 819 " %s "\ 820 % ( e.input(), a.input(), waff), \ 821 ps.oexec) 822 com.run() 823 if (not os.path.isfile(affout) and not ps.dry_run()): 824 self.error_msg("Failed in warping step"); 825 ps.ciao(1) 826 else: 827 self.exists_msg(n.input()) 828 return(affout) 829 830# Main: 831if __name__ == '__main__': 832 833 834 ps = RegWrap('auto_warp.py') 835 ps.init_opts() 836 ps.version() 837 rv = ps.get_user_opts() 838 if (rv != None): ps.ciao(1) 839 840 #process and check input params 841 if(not (ps.process_input())): 842 ps.ciao(1) 843 844 # get rid of any previous temporary data 845 ps.cleanup() 846 847 ps.anat_no_skull = a=ps.prepare_input() 848 b=ps.prepare_base() 849 if (ps.affine_input_xmat == 'AUTO'): #@auto_tlrc 850 if (ps.affine_dxyz != 0.0): 851 a,b = ps.match_resolutions(a,b, '.rs', ps.affine_dxyz) 852 #now run auto_tlrc 853 a, ps.affine_input_xmat = ps.align_auto_tlrc(a,b, at_opt=ps.at_opt) 854 elif (ps.affine_input_xmat == 'ID'): #input already in std space 855 if (0): print("Nothing to do") 856 else: #User specified matrix to take input to std space 857 a, ps.affine_input_xmat = ps.align_auto_tlrc(a,b, xmat=ps.affine_input_xmat) 858 859 #set resolutions properly now with master always the base 860 a,b = ps.match_resolutions(a,b, '.rw', ps.warp_dxyz, m=b) 861 862 #compute the warp 863 a, ps.warp_input_xform = ps.qwarping(a=a, b=b) 864 865 #apply warps 866 # warp datasets may grow, so pass base as master 26 Mar 2014 [rcr/zss] 867 aw = ps.qwarp_applying(a=ps.input, aff=ps.affine_input_xmat, 868 wrp=ps.warp_input_xform, master=b) 869 870 ps.save_history(aw, ps.oexec) 871 872 #apply warps for all followers 873 if (len(ps.followers)): 874 k = 0 875 for fol in ps.followers: 876 if (ps.affine_followers_xmat[k] == 'AUTO'): 877 ps.affine_followers_xmat[k] = ps.align_epi_anat(e=fol, a = ps.anat_no_skull, aff=ps.affine_input_xmat) 878 ps.qwarp_applying(a=fol, aff=ps.affine_followers_xmat[k], wrp=ps.warp_input_xform, dxyz=2.5) 879 k = k + 1 880 881 #cleanup after the parents too? 882 if (ps.rmrm): 883 ps.cleanup() 884 885 ps.ciao(0) 886