antspyt1w.get_data
Get local ANTsPyT1w data
1""" 2Get local ANTsPyT1w data 3""" 4 5__all__ = ['get_data','map_segmentation_to_dataframe','hierarchical', 6 'random_basis_projection', 'deep_dkt','deep_hippo','deep_tissue_segmentation', 7 'deep_brain_parcellation', 'deep_mtl', 'label_hemispheres','brain_extraction', 8 'hemi_reg', 'region_reg', 't1_hypointensity', 'zoom_syn', 'merge_hierarchical_csvs_to_wide_format', 9 'map_intensity_to_dataframe', 'deep_nbm', 'map_cit168', 'deep_cit168','minimal_sr_preprocessing'] 10 11from pathlib import Path 12import os 13from os.path import exists 14import pandas as pd 15import math 16import os.path 17from os import path 18import pickle 19import sys 20import numpy as np 21import random 22import re 23import functools 24from operator import mul 25from scipy.sparse.linalg import svds 26from PyNomaly import loop 27import scipy as sp 28import matplotlib.pyplot as plt 29from PIL import Image 30import scipy.stats as ss 31import warnings 32import ants 33import antspynet 34import tensorflow as tf 35from multiprocessing import Pool 36 37DATA_PATH = os.path.expanduser('~/.antspyt1w/') 38# Internal global variable to hold the seed 39_GLOBAL_SCIENTIFIC_SEED = None 40 41def set_global_scientific_computing_random_seed(seed: int = 42, deterministic: bool = True): 42 """ 43 Set a global seed for reproducibility across Python, NumPy, and scientific ML frameworks. 44 45 Parameters 46 ---------- 47 seed : int 48 The seed to use for all random number generators. 49 deterministic : bool 50 If True, sets flags for deterministic behavior in supported libraries (PyTorch, TensorFlow). 51 """ 52 global _GLOBAL_SCIENTIFIC_SEED 53 _GLOBAL_SCIENTIFIC_SEED = seed 54 55 # Python's built-in RNG 56 random.seed(seed) 57 58 # NumPy RNG 59 np.random.seed(seed) 60 61 # Set Python hash seed for reproducibility 62 os.environ["PYTHONHASHSEED"] = str(seed) 63 64 # TensorFlow 65 try: 66 import tensorflow as tf 67 tf.random.set_seed(seed) 68 if deterministic: 69 os.environ["TF_DETERMINISTIC_OPS"] = "1" 70 tf.config.experimental.enable_op_determinism() 71 except ImportError: 72 pass 73 74 # PyTorch 75 try: 76 import torch 77 torch.manual_seed(seed) 78 torch.cuda.manual_seed(seed) 79 torch.cuda.manual_seed_all(seed) 80 if deterministic: 81 torch.backends.cudnn.deterministic = True 82 torch.backends.cudnn.benchmark = False 83 except ImportError: 84 pass 85 86 # JAX 87 try: 88 import jax 89 jax_key = jax.random.PRNGKey(seed) 90 globals()["_jax_key"] = jax_key 91 except ImportError: 92 pass 93 94 warnings.warn("Remember to set 'random_state=seed' in scikit-learn models.", stacklevel=2) 95 print(f"[INFO] Global scientific computing seed set to {seed}") 96 97 98def get_global_scientific_computing_random_seed(): 99 """ 100 Retrieve the most recently set global seed. 101 102 Returns 103 ------- 104 int or None 105 The seed set by `set_global_scientific_computing_seed()`, or None if not yet set. 106 """ 107 return _GLOBAL_SCIENTIFIC_SEED 108 109set_global_scientific_computing_random_seed(1234) 110 111def get_data(name=None, force_download=False, version=46, target_extension='.csv'): 112 """ 113 Get ANTsPyT1w data filename 114 115 The first time this is called, it will download data to ~/.antspyt1w. 116 After, it will just read data from disk. The ~/.antspyt1w may need to 117 be periodically deleted in order to ensure data is current. 118 119 Arguments 120 --------- 121 name : string 122 name of data tag to retrieve 123 Options: 124 - 'all' 125 - 'dkt' 126 - 'hemisphere' 127 - 'lobes' 128 - 'tissues' 129 - 'T_template0' 130 - 'T_template0_LR' 131 - 'T_template0_LobesBstem' 132 - 'T_template0_WMP' 133 - 'T_template0_Symmetric' 134 - 'T_template0_SymmetricLR' 135 - 'PPMI-3803-20120814-MRI_T1-I340756' 136 - 'simwmseg' 137 - 'simwmdisc' 138 - 'wmh_evidence' 139 - 'wm_major_tracts' 140 141 force_download: boolean 142 143 version: version of data to download (integer) 144 145 Returns 146 ------- 147 string 148 filepath of selected data 149 150 Example 151 ------- 152 >>> import ants 153 >>> ppmi = ants.get_ants_data('ppmi') 154 """ 155 import os 156 import shutil 157 from pathlib import Path 158 import tensorflow as tf 159 DATA_PATH = os.path.join(os.path.expanduser('~'), '.antspyt1w') 160 os.makedirs(DATA_PATH, exist_ok=True) 161 162 def mv_subfolder_files(folder, verbose=False): 163 """ 164 Move files from subfolders to the parent folder. 165 166 Parameters 167 ---------- 168 folder : str 169 Path to the folder. 170 verbose : bool, optional 171 Print information about the files and folders being processed (default is False). 172 173 Returns 174 ------- 175 None 176 """ 177 import os 178 import shutil 179 for root, dirs, files in os.walk(folder): 180 if verbose: 181 print(f"Processing directory: {root}") 182 print(f"Subdirectories: {dirs}") 183 print(f"Files: {files}") 184 185 for file in files: 186 if root != folder: 187 if verbose: 188 print(f"Moving file: {file} from {root} to {folder}") 189 shutil.move(os.path.join(root, file), folder) 190 191 for dir in dirs: 192 if root != folder: 193 if verbose: 194 print(f"Removing directory: {dir} from {root}") 195 shutil.rmtree(os.path.join(root, dir)) 196 197 def download_data(version): 198 url = "https://ndownloader.figshare.com/articles/14766102/versions/" + str(version) 199 target_file_name = "14766102.zip" 200 target_file_name_path = tf.keras.utils.get_file(target_file_name, url, 201 cache_subdir=DATA_PATH, extract=True) 202 os.remove(os.path.join(DATA_PATH, target_file_name)) 203 204 if force_download: 205 download_data(version=version) 206 207 mv_subfolder_files( os.path.expanduser("~/.antspyt1w"), False ) 208 209 # Move files from subdirectories to the main directory 210 for root, dirs, files in os.walk(DATA_PATH): 211 for file in files: 212 if root != DATA_PATH: 213 shutil.move(os.path.join(root, file), DATA_PATH) 214 for dir in dirs: 215 if root != DATA_PATH: 216 shutil.rmtree(os.path.join(root, dir)) 217 218 files = [] 219 for fname in os.listdir(DATA_PATH): 220 if fname.endswith(target_extension): 221 fname = os.path.join(DATA_PATH, fname) 222 files.append(fname) 223 224 if len(files) == 0: 225 download_data(version=version) 226 for fname in os.listdir(DATA_PATH): 227 if fname.endswith(target_extension): 228 fname = os.path.join(DATA_PATH, fname) 229 files.append(fname) 230 231 mv_subfolder_files( os.path.expanduser("~/.antspyt1w"), False ) 232 233 if name == 'all': 234 return files 235 236 datapath = None 237 238 for fname in os.listdir(DATA_PATH): 239 mystem = Path(fname).resolve().stem 240 mystem = Path(mystem).resolve().stem 241 mystem = Path(mystem).resolve().stem 242 if name == mystem and fname.endswith(target_extension): 243 datapath = os.path.join(DATA_PATH, fname) 244 245 if datapath is None: 246 os.listdir(DATA_PATH) 247 warnings.warn("datapath in get_data is None - must be some issue in downloading the data from figshare. try calling get_data in the setup of your system.") 248 return datapath 249 250 251def map_segmentation_to_dataframe( segmentation_type, segmentation_image ): 252 """ 253 Match the segmentation to its appropriate data frame. We do not check 254 if the segmentation_type and segmentation_image match; this may be indicated 255 by the number of missing values on output eg in column VolumeInMillimeters. 256 257 Arguments 258 --------- 259 segmentation_type : string 260 name of segmentation_type data frame to retrieve 261 Options: 262 - 'dkt' 263 - 'lobes' 264 - 'tissues' 265 - 'hemisphere' 266 - 'wm_major_tracts' 267 268 segmentation_image : antsImage with same values (or mostly the same) as are 269 expected by segmentation_type 270 271 Returns 272 ------- 273 dataframe 274 275 """ 276 mydf_fn = get_data( segmentation_type ) 277 mydf = pd.read_csv( mydf_fn ) 278 mylgo = ants.label_geometry_measures( segmentation_image ) 279 return pd.merge( mydf, mylgo, how='left', on=["Label"] ) 280 281 282def map_intensity_to_dataframe(segmentation_type, intensity_image, segmentation_image): 283 """ 284 Match intensity values within segmentation labels to an appropriate data frame. 285 286 Arguments 287 --------- 288 segmentation_type : string or pandas.DataFrame 289 If string: 290 Name of the segmentation type data frame to retrieve using get_data. 291 Options include data available via get_data (e.g., "lobes"). 292 If pandas.DataFrame: 293 A pre-existing DataFrame containing label information. 294 295 intensity_image : ants.ANTsImage 296 ANTsImage with intensity values to summarize. 297 298 segmentation_image : ants.ANTsImage 299 ANTsImage with labels (or mostly the same values) as expected by segmentation_type. 300 301 Returns 302 ------- 303 pandas.DataFrame 304 A DataFrame summarizing the intensity values within segmentation labels. 305 """ 306 if isinstance(segmentation_type, str): 307 # Retrieve the data frame based on the segmentation type string 308 mydf_fn = get_data(segmentation_type) 309 mydf = pd.read_csv(mydf_fn) 310 elif isinstance(segmentation_type, pd.DataFrame): 311 # Use the provided DataFrame 312 mydf = segmentation_type 313 else: 314 raise ValueError("segmentation_type must be either a string or a pandas DataFrame") 315 316 # Get label statistics for the intensity and segmentation images 317 mylgo = ants.label_stats(intensity_image, segmentation_image) 318 319 # Rename the label column to match expected name for merging 320 mylgo = mylgo.rename(columns={'LabelValue': 'Label'}) 321 322 # Merge the segmentation DataFrame with the label statistics 323 result = pd.merge(mydf, mylgo, how='left', on=["Label"]) 324 325 return result 326 327 328def myproduct(lst): 329 return( functools.reduce(mul, lst) ) 330 331def mahalanobis_distance( x ): 332 """ 333 Calculate mahalanobis distance from a dataframe 334 335 Arguments 336 --------- 337 x : dataframe of random projections or other data 338 339 Returns 340 ------- 341 dictionary of distances and outlierness categories 342 343 """ 344 # M-Distance 345 x_minus_mu = x - np.mean(x) 346 cov = np.cov(x.values.T) #Covariance 347 inv_covmat = sp.linalg.inv(cov) #Inverse covariance 348 left_term = np.dot(x_minus_mu, inv_covmat) 349 mahal = np.dot(left_term, x_minus_mu.T) 350 md = np.sqrt(mahal.diagonal()) 351 #Flag as outlier 352 outlier = [] 353 #Cut-off point 354 C = np.sqrt(sp.stats.chi2.ppf((1-0.001), df=x.shape[1])) #degrees of freedom = number of variables 355 for index, value in enumerate(md): 356 if value > C: 357 outlier.append(index) 358 else: 359 continue 360 return { "distance": md, "outlier": outlier } 361 362 363def patch_eigenvalue_ratio( x, n, radii, evdepth = 0.9, mask=None, standardize=False ): 364 """ 365 Patch-based eigenvalue ratio calculation. Computes the ratio (percent) of the eigenvalues needed to reach evdepth explained variance. 366 367 Arguments 368 --------- 369 x : input image - if t1, likely better if brain extracted 370 371 n : number of patches 372 373 radii : list of radius values 374 375 evdepth : value in zero to one 376 377 mask : optional antsImage 378 379 standardize : boolean standardizes the patch matrix if True (subtract mean) 380 381 Returns 382 ------- 383 an eigenvalue ratio in the range of [0,1] 384 385 """ 386 nptch=n 387 radder=radii 388 if mask is None: 389 msk=ants.threshold_image( x, "Otsu", 1 ) 390 else: 391 msk = mask.clone() 392 rnk=ants.rank_intensity(x,msk,True) 393 npatchvox = myproduct( radder ) 394 ptch0 = ants.extract_image_patches( rnk, tuple(radder), mask_image=msk, 395 max_number_of_patches = nptch, return_as_array=False, randomize=False ) 396 ptch = [] 397 for k in range(len(ptch0)): 398 if np.prod( ptch0[k].shape ) == npatchvox : 399 ptch.append( np.reshape( ptch0[k], npatchvox ) ) 400 X = np.stack( ptch ) 401 if standardize: 402 X = X - X.mean(axis=0, keepdims=True) 403 # u, s, v = svds(X , min(X.shape)-1 ) 404 thespectrum = np.linalg.svd( X, compute_uv=False ) 405 spectralsum = thespectrum.sum() 406 targetspec = spectralsum * evdepth 407 spectralcumsum = np.cumsum( thespectrum ) 408 numer = np.argmin( abs( spectralcumsum - evdepth * spectralsum ) ) 409 denom = len( thespectrum ) 410 # print( str(numer) + " " + str(denom) ) 411 return numer/denom 412 413def loop_outlierness( random_projections, reference_projections=None, 414 standardize=True, extent=3, n_neighbors=24, cluster_labels=None ): 415 """ 416 Estimate loop outlierness for the input. 417 418 Arguments 419 --------- 420 random_projections : output of random_basis_projection or a similar function 421 422 reference_projections : produced in the same way as random_projections 423 424 standardize: boolean will subtract mean and divide by sd of the reference data 425 426 extent: an integer value [1, 2, 3] that controls the statistical 427 extent, e.g. lambda times the standard deviation from the mean (optional, 428 default 2) 429 430 n_neighbors: the total number of neighbors to consider w.r.t. each 431 sample (optional, default 10) 432 433 cluster_labels: a numpy array of cluster assignments w.r.t. each 434 sample (optional, default None) 435 436 Returns 437 ------- 438 loop outlierness probability 439 440 """ 441 use_reference=True 442 if reference_projections is None: 443 use_reference = False 444 reference_projections = random_projections 445 nBasisUse = reference_projections.shape[1] 446 if random_projections.shape[1] < nBasisUse: 447 nBasisUse = random_projections.shape[1] 448 449 refbases = reference_projections.iloc[:,:nBasisUse] 450 myax=0 451 refbasesmean = refbases.mean(axis=myax) 452 refbasessd = refbases.std(axis=myax) 453 normalized_df = refbases 454 if standardize: 455 normalized_df = (normalized_df-refbasesmean)/refbasessd 456 if use_reference: 457 temp = random_projections.iloc[:,:nBasisUse] 458 if standardize: 459 temp = (temp-refbasesmean)/refbasessd 460 normalized_df = pd.concat( [normalized_df, temp], axis=0 ).dropna(axis=0) 461 if cluster_labels is None: 462 m = loop.LocalOutlierProbability(normalized_df, 463 extent=extent, 464 n_neighbors=n_neighbors ).fit() 465 else: 466 m = loop.LocalOutlierProbability(normalized_df, 467 extent=extent, 468 n_neighbors=n_neighbors, 469 cluster_labels=cluster_labels).fit() 470 scores = m.local_outlier_probabilities 471 472 return scores 473 474 475 476def random_basis_projection( x, template, 477 type_of_transform='Similarity', 478 refbases = None, 479 nBasis=10, random_state = 99 ): 480 """ 481 Produce unbiased data descriptors for a given image which can be used 482 to assist data inspection and ranking. can be used with any image 483 brain extracted or not, any modality etc. but we assume we can 484 meaningfully map to a template, at least with a low-dimensional 485 transformation, e.g. Translation, Rigid, Similarity. 486 487 Arguments 488 --------- 489 x : antsImage 490 491 template : antsImage reference template 492 493 type_of_transform: one of Translation, Rigid, Similarity, Affine 494 495 refbases : reference bases for outlierness calculations 496 497 nBasis : number of variables to derive 498 499 random_state : seed 500 501 Returns 502 ------- 503 dataframe with projections and an outlierness estimate. 504 505 the outlierness estimate is based on a small reference dataset of young controls. 506 the user/researcher may want to use a different reference set. see the 507 function loop_outlierness for one way to do that. 508 509 """ 510 template = ants.crop_image( template ) 511 template = ants.iMath( template, "Normalize" ) 512 np.random.seed(int(random_state)) 513 nvox = template.shape 514 # X = np.random.rand( nBasis+1, myproduct( nvox ) ) 515 # u, s, randbasis = svds(X, k=nBasis) 516 # if randbasis.shape[1] != myproduct(nvox): 517 # raise ValueError("columns in rand basis do not match the nvox product") 518 519 randbasis = np.random.randn( myproduct( nvox ), nBasis ) 520 rbpos = randbasis.copy() 521 rbpos[rbpos<0] = 0 522 norm = ants.iMath( x, "Normalize" ) 523 trans = ants.registration( template, norm, 524 type_of_transform='antsRegistrationSyNQuickRepro[t]' ) 525 resamp = ants.registration( template, norm, 526 type_of_transform=type_of_transform, total_sigma=0.5, 527 # aff_metric='GC', 528 random_seed=1, initial_transform=trans['fwdtransforms'][0] )['warpedmovout'] 529# mydelta = ants.from_numpy( ( resamp - template ).abs() ) 530 mydelta = resamp - template 531 imat = ants.get_neighborhood_in_mask( mydelta, mydelta*0+1,[0,0,0], boundary_condition='mean' ) 532 uproj = np.matmul(imat, randbasis) 533 uprojpos = np.matmul(imat, rbpos) 534 record = {} 535 uproj_counter = 0 536 for i in uproj[0]: 537 uproj_counter += 1 538 name = "RandBasisProj" + str(uproj_counter).zfill(2) 539 record[name] = i 540 uprojpos_counter = 0 541 for i in uprojpos[0]: 542 uprojpos_counter += 1 543 name = "RandBasisProjPos" + str(uprojpos_counter).zfill(2) 544 record[name] = i 545 df = pd.DataFrame(record, index=[0]) 546 547 if refbases is None: 548 refbases = pd.read_csv( get_data( "reference_basis", target_extension='.csv' ) ) 549 df['loop_outlier_probability'] = loop_outlierness( df, refbases, 550 n_neighbors=refbases.shape[0] )[ refbases.shape[0] ] 551 mhdist = 0.0 552 if nBasis == 10: 553 temp = pd.concat( [ refbases, df.iloc[:,:nBasis] ], axis=0 ) 554 mhdist = mahalanobis_distance( temp )['distance'][ refbases.shape[0] ] 555 df['mhdist'] = mhdist 556 df['templateL1']=mydelta.abs().mean() 557 return df 558 559 560def resnet_grader( x, weights_filename = None ): 561 """ 562 Supervised grader / scoring of t1 brain 563 564 Arguments 565 --------- 566 567 x : antsImage of t1 brain 568 569 weights_filename : optional weights filename 570 571 Returns 572 ------- 573 two letter grades 574 575 """ 576 577 if weights_filename is None: 578 weights_filename=get_data( 'resnet_grader', target_extension='.h5' ) 579 580 if not exists( weights_filename ): 581 print("resnet_grader weights do not exist: " + weights_filename ) 582 return None 583 584 mdl = antspynet.create_resnet_model_3d( [None,None,None,1], 585 lowest_resolution = 32, 586 number_of_outputs = 4, 587 cardinality = 1, 588 squeeze_and_excite = False ) 589 mdl.load_weights( weights_filename ) 590 591 592 t1 = ants.iMath( x - x.min(), "Normalize" ) 593 bxt = ants.threshold_image( t1, 0.01, 1.0 ) 594 t1 = ants.rank_intensity( t1, mask=bxt, get_mask=True ) 595 templateb = ants.image_read( get_data( "S_template3_brain", target_extension='.nii.gz' ) ) 596 templateb = ants.crop_image( templateb ).resample_image( [1,1,1] ) 597 templateb = ants.pad_image_by_factor( templateb, 8 ) 598 templatebsmall = ants.resample_image( templateb, [2,2,2] ) 599 reg = ants.registration( templatebsmall, t1, 'Similarity', verbose=False ) 600 ilist = list() 601 refimg=templateb 602 ilist.append( [refimg] ) 603 nsim = 16 604 uu = ants.randomly_transform_image_data( refimg, ilist, 605 number_of_simulations = nsim, 606 transform_type='scaleShear', sd_affine=0.075 ) 607 fwdaffgd = ants.read_transform( reg['fwdtransforms'][0]) 608 scoreNums = np.zeros( 4 ) 609 scoreNums[3]=0 610 scoreNums[2]=1 611 scoreNums[1]=2 612 scoreNums[0]=3 613 scoreNums=scoreNums.reshape( [4,1] ) 614 615 def get_grade( score, probs ): 616 grade='f' 617 if score >= 2.25: 618 grade='a' 619 elif score >= 1.5: 620 grade='b' 621 elif score >= 0.75: 622 grade='c' 623 probgradeindex = np.argmax( probs ) 624 probgrade = ['a','b','c','f'][probgradeindex] 625 return [grade, probgrade, float( score )] 626 627 gradelistNum = [] 628 gradelistProb = [] 629 gradeScore = [] 630 for k in range( nsim ): 631 simtx = uu['simulated_transforms'][k] 632 cmptx = ants.compose_ants_transforms( [simtx,fwdaffgd] ) # good 633 subjectsim = ants.apply_ants_transform_to_image( cmptx, t1, refimg, interpolation='linear' ) 634 subjectsim = ants.add_noise_to_image( subjectsim, 'additivegaussian', (0,0.01) ) 635 xarr = subjectsim.numpy() 636 newshape = list( xarr.shape ) 637 newshape = [1] + newshape + [1] 638 xarr = np.reshape( xarr, newshape ) 639 preds = mdl.predict( xarr ) 640 predsnum = tf.matmul( preds, scoreNums ) 641 locgrades = get_grade( predsnum, preds ) 642 gradelistNum.append( locgrades[0] ) 643 gradelistProb.append( locgrades[1] ) 644 gradeScore.append( locgrades[2] ) 645 646 def most_frequent(List): 647 return max(set(List), key = List.count) 648 649 mydf = pd.DataFrame( { 650 "NumericGrade": gradelistNum, 651 "ProbGrade": gradelistProb, 652 "NumericScore": gradeScore, 653 'grade': most_frequent( gradelistProb ) 654 }) 655 656 smalldf = pd.DataFrame( { 657 'gradeLetter': [mydf.grade[0]], 658 'gradeNum': [mydf.NumericScore.mean()] 659 }, index=[0] ) 660 # print( mydf.Num.value_counts() ) 661 # print( mydf.Prob.value_counts() ) 662 return smalldf 663 664 665 666def inspect_raw_t1( x, output_prefix, option='both' ): 667 """ 668 Quick inspection and visualization of a raw T1 whole head image, 669 whole head and brain or both. The reference data was developed with 670 option both and this will impact results. For the head image, outlierness 671 is estimated vi a quick extraction from background. Variability in this 672 extraction (Otsu) may lead to reliability issues. 673 674 Arguments 675 --------- 676 677 x : antsImage of t1 whole head 678 679 output_prefix: a path and prefix for outputs 680 681 option : string both, brain or head 682 683 Returns 684 ------- 685 two dataframes (head, brain) with projections and outlierness estimates. 686 687 """ 688 689 if x.dimension != 3: 690 raise ValueError('inspect_raw_t1: input image should be 3-dimensional') 691 692 x = ants.iMath( x, "Normalize" ) 693 csvfn = output_prefix + "head.csv" 694 pngfn = output_prefix + "head.png" 695 csvfnb = output_prefix + "brain.csv" 696 pngfnb = output_prefix + "brain.png" 697 698 # reference bases 699 rbh = pd.read_csv( get_data( "refbasis_head", target_extension=".csv" ) ) 700 rbb = pd.read_csv( get_data( "refbasis_brain", target_extension=".csv" ) ) 701 702 # whole head outlierness 703 rbp=None 704 if option == 'both' or option == 'head': 705 bfn = antspynet.get_antsxnet_data( "S_template3" ) 706 templateb = ants.image_read( bfn ).iMath("Normalize") 707 templatesmall = ants.resample_image( templateb, (2,2,2), use_voxels=False ) 708 lomask = ants.threshold_image( x, "Otsu", 2 ).threshold_image(1,2) 709 t1 = ants.rank_intensity( x * lomask, mask=lomask, get_mask=False ) 710 ants.plot( t1, axis=2, nslices=21, ncol=7, filename=pngfn, crop=True ) 711 rbp = random_basis_projection( t1, templatesmall, 712 type_of_transform='Rigid', 713 refbases=rbh ) 714 rbp.to_csv( csvfn ) 715 # fix up the figure 716 looper=float(rbp['loop_outlier_probability'].iloc[0]) 717 ttl="LOOP: " + "{:0.4f}".format(looper) + " MD: " + "{:0.4f}".format(float(rbp['mhdist'].iloc[0])) 718 img = Image.open( pngfn ).copy() 719 plt.figure(dpi=300) 720 plt.imshow(img) 721 plt.text(20, 0, ttl, color="red", fontsize=12 ) 722 plt.axis("off") 723 plt.subplots_adjust(0,0,1,1) 724 plt.savefig( pngfn, bbox_inches='tight',pad_inches = 0) 725 plt.close() 726 727 # same for brain 728 rbpb=None 729 evratio=None 730 if option == 'both' or option == 'brain': 731 if option == 'both': 732 t1 = ants.iMath( x, "TruncateIntensity",0.001, 0.999).iMath("Normalize") 733 lomask = antspynet.brain_extraction( t1, "t1" ) 734 t1 = ants.rank_intensity( t1 * lomask, get_mask=True ) 735 else: 736 t1 = ants.iMath( x, "Normalize" ) 737 t1 = ants.rank_intensity( t1, get_mask=True ) 738 ants.plot( t1, axis=2, nslices=21, ncol=7, filename=pngfnb, crop=True ) 739 templateb = ants.image_read( get_data( "S_template3_brain", target_extension='.nii.gz' ) ) 740 templatesmall = ants.resample_image( templateb, (2,2,2), use_voxels=False ) 741 rbpb = random_basis_projection( t1, 742 templatesmall, 743 type_of_transform='Rigid', 744 refbases=rbb ) 745 rbpb['evratio'] = patch_eigenvalue_ratio( t1, 512, [20,20,20], evdepth = 0.9 ) 746 grade0 = resnet_grader( t1 ).gradeNum[0] 747 msk=ants.threshold_image(t1,0.01,1.0) 748 t1tx=ants.n4_bias_field_correction( t1, mask=msk ) 749 t1tx=ants.iMath(t1tx,'TruncateIntensity',0.001,0.98) 750 grade1 = resnet_grader( t1tx ).gradeNum[0] 751 if grade1 > grade0: 752 grade0 = grade1 753 rbpb['resnetGrade'] = grade0 754 rbpb.to_csv( csvfnb ) 755 looper = float(rbpb['loop_outlier_probability'].iloc[0]) 756 myevr = float( rbpb['evratio'].iloc[0] ) 757 mygrd = float( rbpb['resnetGrade'].iloc[0] ) 758 myl1 = float( rbpb['templateL1'].iloc[0] ) 759 ttl="LOOP: " + "{:0.4f}".format(looper) + " MD: " + "{:0.4f}".format(float(rbpb['mhdist'].iloc[0])) + " EVR: " + "{:0.4f}".format(myevr) + " TL1: " + "{:0.4f}".format(myl1) + " grade: " + "{:0.4f}".format(mygrd) 760 img = Image.open( pngfnb ).copy() 761 plt.figure(dpi=300) 762 plt.imshow(img) 763 plt.text(20, 0, ttl, color="red", fontsize=12 ) 764 plt.axis("off") 765 plt.subplots_adjust(0,0,1,1) 766 plt.savefig( pngfnb, bbox_inches='tight',pad_inches = 0) 767 plt.close() 768 769 return { 770 "head": rbp, 771 "head_image": pngfn, 772 "brain": rbpb, 773 "brain_image": pngfnb, 774 } 775 776def icv( x ): 777 """ 778 Intracranial volume (ICV) estimation from a brain extracted image 779 780 Arguments 781 --------- 782 x : antsImage raw t1 783 784 Returns 785 ------- 786 ICV in cubic millimeters 787 788 """ 789 icvseg = antspynet.brain_extraction( ants.iMath( x, "Normalize" ), modality="t1threetissue")['segmentation_image'].threshold_image(1,2) 790 np.product = np.prod( ants.get_spacing( x ) ) * icvseg.sum() 791 return pd.DataFrame({'icv': [np.product]}) 792 793def brain_extraction( x, dilation = 8.0, method = 'v1', deform=True, verbose=False ): 794 """ 795 quick brain extraction for individual images 796 797 x: input image 798 799 dilation: amount to dilate first brain extraction in millimeters 800 801 method: version currently v0 or any other string gives two different results 802 803 deform: map the image to the training hypersphere before running the extraction 804 805 verbose: boolean 806 807 """ 808 return antspynet.brain_extraction( ants.iMath( x, "Normalize" ), modality="t1threetissue")['segmentation_image'].threshold_image(1,1) 809 810 if deform: 811 reorient_template_file_name_path = antspynet. get_antsxnet_data("S_template3" ) 812 template = ants.image_read( reorient_template_file_name_path ) 813 reg = ants.registration( template, x, 'antsRegistrationSyNQuickRepro[s]', 814 total_sigma=0.5 ) 815 816 closedilmm = 5.0 817 spacing = ants.get_spacing(x) 818 spacing_product = spacing[0] * spacing[1] * spacing[2] 819 spcmin = min( spacing ) 820 dilationRound = int(np.round( dilation / spcmin )) 821 closedilRound = int(np.round( closedilmm / spcmin )) 822 if deform: 823 xn3 = ants.n3_bias_field_correction( reg['warpedmovout'], 8 ).n3_bias_field_correction( 4 ) 824 xn3 = ants.iMath(xn3, "TruncateIntensity",0.001,0.999).iMath("Normalize") 825 else: 826 xn3 = ants.n3_bias_field_correction( x, 8 ).n3_bias_field_correction( 4 ) 827 xn3 = ants.iMath(xn3, "TruncateIntensity",0.001,0.999).iMath("Normalize") 828 if method == 'v0': 829 if verbose: 830 print("method v0") 831 bxtmethod = 't1combined[' + str(closedilRound) +']' # better for individual subjects 832 bxt = antspynet.brain_extraction( xn3, bxtmethod ).threshold_image(2,3).iMath("GetLargestComponent",50).iMath("FillHoles") 833 if deform: 834 bxt = ants.apply_transforms( x, bxt, reg['invtransforms'], interpolator='nearestNeighbor' ) 835 return bxt 836 if method == 'v1': 837 if verbose: 838 print("method v1") 839 bxt = antspynet.brain_extraction( xn3, 't1' ).threshold_image(0.5,3).iMath("GetLargestComponent",50).iMath("FillHoles") 840 if deform: 841 bxt = ants.apply_transforms( x, bxt, reg['invtransforms'], interpolator='nearestNeighbor' ) 842 return bxt 843 if verbose: 844 print("method candidate") 845 bxt0 = antspynet.brain_extraction( xn3, "t1" ).threshold_image(0.5,1.0).iMath("GetLargestComponent",50).morphology( "close", closedilRound ).iMath("FillHoles") 846 bxt0dil = ants.iMath( bxt0, "MD", dilationRound ) 847 image = ants.iMath( xn3 * bxt0dil,"Normalize")*255 848 # no no brainer 849 model = antspynet.create_nobrainer_unet_model_3d((None, None, None, 1)) 850 weights_file_name = antspynet.get_pretrained_network("brainExtractionNoBrainer") 851 model.load_weights(weights_file_name) 852 image_array = image.numpy() 853 image_robust_range = np.quantile(image_array[np.where(image_array != 0)], (0.02, 0.98)) 854 threshold_value = 0.10 * (image_robust_range[1] - image_robust_range[0]) + image_robust_range[0] 855 thresholded_mask = ants.threshold_image(image, -10000, threshold_value, 0, 1) 856 thresholded_image = image * thresholded_mask 857 image_resampled = ants.resample_image(thresholded_image, (256, 256, 256), use_voxels=True) 858 image_array = np.expand_dims(image_resampled.numpy(), axis=0) 859 image_array = np.expand_dims(image_array, axis=-1) 860 brain_mask_array = np.squeeze(model.predict(image_array)) 861 brain_mask_resampled = ants.copy_image_info(image_resampled, ants.from_numpy(brain_mask_array)) 862 brain_mask_image = ants.resample_image(brain_mask_resampled, image.shape, use_voxels=True, interp_type=1) 863 brain_mask_image = brain_mask_image * ants.threshold_image( brain_mask_image, 0.5, 1e9 ) 864 minimum_brain_volume = round(649933.7/spacing_product) 865 bxt1 = ants.label_clusters(brain_mask_image, minimum_brain_volume) 866 nblabels = np.unique( bxt1.numpy() )# .astype(int) 867 maxol = 0.0 868 bestlab = bxt0 869 870 for j in range(1,len(nblabels)): 871 temp = ants.threshold_image( bxt1, j, j ) 872 tempsum = ( bxt0 * temp ).sum() 873 dicer = ants.label_overlap_measures( temp, bxt0 ) 874 if tempsum > maxol and dicer['MeanOverlap'][1] > 0.5 : 875 if verbose: 876 print( str(j) + ' dice ' + str( dicer['MeanOverlap'][1] ) ) 877 maxol = tempsum 878 bestlab = temp 879 adder = trim_segmentation_by_distance( bxt0, 1, closedilmm ) 880 bestlab = ants.threshold_image( bestlab + adder, 1, 2 ) 881 return bestlab 882 883 884 885def subdivide_labels( x, verbose = False ): 886 """ 887 quick subdivision of the labels in a label image. can be applied recursively 888 to get several levels of subdivision. 889 890 x: input hemisphere label image from label_hemispheres 891 892 verbose: boolean 893 894 Example: 895 896 > x=ants.image_read( ants.get_data("ch2") ) 897 > seg=ants.threshold_image(x,1,np.math.inf) 898 > for n in range(5): 899 > seg=subdivide_labels(seg,verbose=True) 900 > ants.plot(x,seg,axis=2,nslices=21,ncol=7,crop=True) 901 902 """ 903 notzero=ants.threshold_image( x, 1, np.math.inf ) 904 ulabs = np.unique( x.numpy() ) 905 ulabs.sort() 906 newx = x * 0.0 907 for u in ulabs: 908 if u > 0: 909 temp = ants.threshold_image( x, u, u ) 910 subimg = ants.crop_image( x, temp ) * 0 911 localshape=subimg.shape 912 axtosplit=np.argmax(localshape) 913 mid=int(np.round( localshape[axtosplit] /2 )) 914 nextlab = newx.max()+1 915 if verbose: 916 print( "label: " + str( u ) ) 917 print( subimg ) 918 print( nextlab ) 919 if axtosplit == 1: 920 subimg[:,0:mid,:]=subimg[:,0:mid,:]+nextlab 921 subimg[:,(mid):(localshape[axtosplit]),:]=subimg[:,(mid):(localshape[axtosplit]),:]+nextlab+1 922 if axtosplit == 0: 923 subimg[0:mid,:,:]=subimg[0:mid,:,:]+nextlab 924 subimg[(mid):(localshape[axtosplit]),:,:]=subimg[(mid):(localshape[axtosplit]),:,:]+nextlab+1 925 if axtosplit == 2: 926 subimg[:,:,0:mid]=subimg[:,:,0:mid]+nextlab 927 subimg[:,:,(mid):(localshape[axtosplit])]=subimg[:,:,(mid):(localshape[axtosplit])]+nextlab+1 928 newx = newx + ants.resample_image_to_target( subimg, newx, interp_type='nearestNeighbor' ) * notzero 929 return newx 930 931 932 933 934def subdivide_hemi_label( x ): 935 """ 936 quick subdivision of the hemisphere label to go from a 2-label to a 4-label. 937 will subdivide along the largest axis in terms of voxels. 938 939 x: input hemisphere label image from label_hemispheres 940 941 verbose: boolean 942 943 """ 944 notzero=ants.threshold_image( x, 1, 1e9 ) 945 localshape=ants.crop_image( x, ants.threshold_image( x, 1, 1 ) ).shape 946 axtosplit=np.argmax(localshape) 947 mid=int(np.round( localshape[axtosplit] /2 )) 948 if axtosplit == 1: 949 x[:,0:mid,:]=x[:,0:mid,:]+3 950 x[:,(mid):(localshape[axtosplit]),:]=x[:,(mid):(localshape[axtosplit]),:]+5 951 if axtosplit == 0: 952 x[0:mid,:,:]=x[0:mid,:,:]+3 953 x[(mid):(localshape[axtosplit]),:,:]=x[(mid):(localshape[axtosplit]),:,:]+5 954 if axtosplit == 2: 955 x[:,:,0:mid]=x[:,:,0:mid]+3 956 x[:,:,(mid):(localshape[axtosplit])]=x[:,:,(mid):(localshape[axtosplit])]+5 957 return x*notzero 958 959 960def special_crop( x, pt, domainer ): 961 """ 962 quick cropping to a fixed size patch around a center point 963 964 x: input image 965 966 pt: list of physical space coordinates 967 968 domainer: list defining patch size 969 970 verbose: boolean 971 972 """ 973 pti = np.round( ants.transform_physical_point_to_index( x, pt ) ) 974 xdim = x.shape 975 for k in range(len(xdim)): 976 if pti[k] < 0: 977 pti[k]=0 978 if pti[k] > (xdim[k]-1): 979 pti[k]=(xdim[k]-1) 980 mim = ants.make_image( domainer ) 981 ptioff = pti.copy() 982 for k in range(len(xdim)): 983 ptioff[k] = ptioff[k] - np.round( domainer[k] / 2 ) 984 domainerlo = [] 985 domainerhi = [] 986 for k in range(len(xdim)): 987 domainerlo.append( int(ptioff[k] - 1) ) 988 domainerhi.append( int(ptioff[k] + 1) ) 989 loi = ants.crop_indices( x, tuple(domainerlo), tuple(domainerhi) ) 990 mim = ants.copy_image_info(loi,mim) 991 return ants.resample_image_to_target( x, mim ) 992 993def label_hemispheres( x, template=None, templateLR=None, reg_iterations=[200,50,2,0] ): 994 """ 995 quick somewhat noisy registration solution to hemisphere labeling. typically 996 we label left as 1 and right as 2. 997 998 x: input image 999 1000 template: MNI space template, should be "croppedMni152" or "biobank" 1001 1002 templateLR: a segmentation image of template hemispheres 1003 1004 reg_iterations: reg_iterations for ants.registration 1005 1006 """ 1007 1008 if template is None or templateLR is None: 1009 tfn = get_data('T_template0', target_extension='.nii.gz' ) 1010 tfnw = get_data('T_template0_WMP', target_extension='.nii.gz' ) 1011 tlrfn = get_data('T_template0_LR', target_extension='.nii.gz' ) 1012 bfn = antspynet.get_antsxnet_data( "croppedMni152" ) 1013 template = ants.image_read( tfn ) 1014 template = ( template * antspynet.brain_extraction( template, 't1' ) ).iMath( "Normalize" ) 1015 templateawmprior = ants.image_read( tfnw ) 1016 templateLR = ants.image_read( tlrfn ) 1017 1018 reg = ants.registration( 1019 ants.rank_intensity(x), 1020 ants.rank_intensity(template), 1021 'antsRegistrationSyNQuickRepro[s]', 1022 aff_metric='GC', 1023 syn_metric='CC', 1024 syn_sampling=2, 1025 reg_iterations=reg_iterations, total_sigma=0.5, 1026 random_seed = 1 ) 1027 return( ants.apply_transforms( x, templateLR, reg['fwdtransforms'], 1028 interpolator='nearestNeighbor') ) 1029 1030def deep_tissue_segmentation( x, template=None, registration_map=None, 1031 atropos_prior=None, sr_model=None ): 1032 """ 1033 modified slightly more efficient deep atropos that also handles the 1034 extra CSF issue. returns segmentation and probability images. see 1035 the tissues csv available from get_data. 1036 1037 x: input image, raw 1038 1039 template: MNI space template, should be "croppedMni152" or "biobank" 1040 1041 registration_map: pre-existing output from ants.registration 1042 1043 atropos_prior: prior weight for atropos post-processing; set to None if you 1044 do not want to use this. will modify the CSF, GM and WM segmentation to 1045 better fit the image intensity at the resolution of the input image. 1046 1047 sr_model: optional (FIXME) 1048 1049 """ 1050 1051 dapper = antspynet.deep_atropos( x, do_denoising=False ) 1052 1053 myk='segmentation_image' 1054 if atropos_prior is not None: 1055 msk = ants.threshold_image( dapper['segmentation_image'], 2, 3 ).iMath("GetLargestComponent",50) 1056 msk = ants.morphology( msk, "close", 2 ) 1057 mskfull = ants.threshold_image( dapper['segmentation_image'], 1, 6 ) 1058 mskfull = mskfull - msk 1059 priors = dapper[myk][1:4] 1060 for k in range( len( priors ) ): 1061 priors[k]=ants.image_clone( priors[k]*msk ) 1062 aap = ants.atropos( x, msk, i=priors, m='[0.0,1x1x1]', 1063 c = '[1,0]', priorweight=atropos_prior, verbose=1 ) 1064 dapper['segmentation_image'] = aap['segmentation'] * msk + dapper['segmentation_image'] * mskfull 1065 dapper['probability_images'][1:4] = aap['probabilityimages'] 1066 1067 return dapper 1068 1069def deep_brain_parcellation( 1070 target_image, 1071 template, 1072 img6seg = None, 1073 do_cortical_propagation=False, 1074 atropos_prior=None, 1075 verbose=True, 1076): 1077 """ 1078 modified slightly more efficient deep dkt that also returns atropos output 1079 thus providing a complete hierarchical parcellation of t1w. we run atropos 1080 here so we dont need to redo registration separately. see 1081 the lobes and dkt csv available from get_data. 1082 1083 target_image: input image 1084 1085 template: MNI space template, should be "croppedMni152" or "biobank" 1086 1087 img6seg: optional pre-existing brain segmentation - a 6-class image ( ANTs standard ) 1088 1089 do_cortical_propagation: boolean, adds a bit extra time to propagate cortical 1090 labels explicitly into cortical segmentation 1091 1092 atropos_prior: prior weight for atropos post-processing; set to None if you 1093 do not want to use this. will modify the CSF, GM and WM segmentation to 1094 better fit the image intensity at the resolution of the input image. 1095 1096 verbose: boolean 1097 1098 1099 Returns 1100 ------- 1101 a dictionary containing: 1102 1103 - tissue_segmentation : 6 tissue segmentation 1104 - tissue_probabilities : probability images associated with above 1105 - dkt_parcellation : tissue agnostic DKT parcellation 1106 - dkt_lobes : major lobes of the brain 1107 - dkt_cortex: cortical tissue DKT parcellation (if requested) 1108 - hemisphere_labels: free to get hemisphere labels 1109 - wmSNR : white matter signal-to-noise ratio 1110 - wmcsfSNR : white matter to csf signal-to-noise ratio 1111 1112 """ 1113 if verbose: 1114 print("Begin registration") 1115 1116 rig = ants.registration( template, ants.rank_intensity(target_image), 1117 "antsRegistrationSyNQuickRepro[a]", 1118 aff_iterations = (500,200,0,0), 1119 total_sigma=0.5, 1120 random_seed=1 ) 1121 rigi = ants.apply_transforms( template, target_image, rig['fwdtransforms']) 1122 1123 if verbose: 1124 print("Begin Atropos tissue segmentation") 1125 1126 if img6seg is None: 1127 if verbose: 1128 print( "do deep_tissue_segmentation") 1129 mydap = deep_tissue_segmentation( target_image ) 1130 else: 1131 if verbose: 1132 print( "use existing deep_tissue_segmentation") 1133 mydap = { 'segmentation_image': img6seg, 'probability_images': None } 1134 1135 if verbose: 1136 print("End Atropos tissue segmentation") 1137 1138 if verbose: 1139 print("Begin DKT") 1140 1141 dktprepro = True 1142 dkt = antspynet.desikan_killiany_tourville_labeling( 1143 target_image, 1144 do_preprocessing=True, 1145 return_probability_images=False, 1146 do_lobar_parcellation = True, 1147 do_denoising=False 1148 ) 1149 1150 myhemiL = ants.threshold_image( dkt['lobar_parcellation'], 1, 6 ) 1151 myhemiR = ants.threshold_image( dkt['lobar_parcellation'], 7, 12 ) 1152 myhemi = myhemiL + myhemiR * 2.0 1153 brainmask = ants.threshold_image( mydap['segmentation_image'], 1, 6 ) 1154 myhemi = ants.iMath( brainmask, 'PropagateLabelsThroughMask', myhemi, 100, 0) 1155 1156 cortprop = None 1157 if do_cortical_propagation: 1158 cortprop = ants.threshold_image( mydap['segmentation_image'], 2, 2 ) 1159 cortlab = dkt['segmentation_image'] * ants.threshold_image( dkt['segmentation_image'], 1000, 5000 ) 1160 cortprop = ants.iMath( cortprop, 'PropagateLabelsThroughMask', 1161 cortlab, 1, 0) 1162 1163 wmseg = ants.threshold_image( mydap['segmentation_image'], 3, 3 ) 1164 wmMean = target_image[ wmseg == 1 ].mean() 1165 wmStd = target_image[ wmseg == 1 ].std() 1166 csfseg = ants.threshold_image( mydap['segmentation_image'], 1, 1 ) 1167 csfStd = target_image[ csfseg == 1 ].std() 1168 wmSNR = wmMean/wmStd 1169 wmcsfSNR = wmMean/csfStd 1170 1171 return { 1172 "tissue_segmentation":mydap['segmentation_image'], 1173 "tissue_probabilities":mydap['probability_images'], 1174 "dkt_parcellation":dkt['segmentation_image'], 1175 "dkt_lobes":dkt['lobar_parcellation'], 1176 "dkt_cortex": cortprop, 1177 "hemisphere_labels": myhemi, 1178 "wmSNR": wmSNR, 1179 "wmcsfSNR": wmcsfSNR, } 1180 1181 1182def deep_hippo( 1183 img, 1184 template, 1185 number_of_tries = 10, 1186 tx_type="antsRegistrationSyNQuickRepro[a]", 1187 sr_model = None, 1188 verbose=False 1189): 1190 1191 1192 super_resolution = None 1193 1194 avgleft = img * 0 1195 avgright = img * 0 1196 for k in range(number_of_tries): 1197 if verbose: 1198 print( "try " + str(k)) 1199 rig = ants.registration( template, ants.rank_intensity(img), 1200 tx_type, random_seed=k, total_sigma=0.5, verbose=verbose ) 1201 if verbose: 1202 print( rig ) 1203 rigi = ants.apply_transforms( template, img, rig['fwdtransforms'] ) 1204 if verbose: 1205 print( "done with warp - do hippmapp3r" ) 1206 hipp = antspynet.hippmapp3r_segmentation( rigi, do_preprocessing=False ) 1207 if verbose: 1208 print( "done with hippmapp3r - backtransform" ) 1209 1210 if sr_model is not None: 1211 mysr = super_resolution_segmentation_per_label( 1212 rigi, segmentation=hipp, upFactor=[2,2,2], 1213 sr_model=sr_model, segmentation_numbers=[1,2], 1214 dilation_amount=0, 1215 probability_images=None, 1216 probability_labels=[1,2], max_lab_plus_one=True, verbose=True 1217 ) 1218 hipp = mysr['super_resolution_segmentation'] 1219 1220 hippr = ants.apply_transforms( 1221 img, 1222 hipp, 1223 rig['fwdtransforms'], 1224 whichtoinvert=[True], 1225 interpolator='nearestNeighbor', 1226 ) 1227 avgleft = avgleft + ants.threshold_image( hippr, 2, 2 ).iMath("GetLargestComponent",1) / float(number_of_tries) 1228 avgright = avgright + ants.threshold_image( hippr, 1, 1 ).iMath("GetLargestComponent",1) / float(number_of_tries) 1229 1230 1231 avgright = ants.iMath(avgright,"Normalize") # output: probability image right 1232 avgleft = ants.iMath(avgleft,"Normalize") # output: probability image left 1233 hippright_bin = ants.threshold_image( avgright, 0.5, 2.0 ).iMath("GetLargestComponent",1) 1234 hippleft_bin = ants.threshold_image( avgleft, 0.5, 2.0 ).iMath("GetLargestComponent",1) 1235 hipp_bin = hippleft_bin + hippright_bin * 2 1236 1237 hippleftORlabels = ants.label_geometry_measures(hippleft_bin, avgleft) 1238 hippleftORlabels['Description'] = 'left hippocampus' 1239 hipprightORlabels = ants.label_geometry_measures(hippright_bin, avgright) 1240 hipprightORlabels['Description'] = 'right hippocampus' 1241 # hippleftORlabels=hippleftORlabels.append( hipprightORlabels ) 1242 # FIXME possible fix below 1243 hippleftORlabels = pd.concat( [hippleftORlabels, hipprightORlabels] , axis=1 ) 1244 hippleftORlabels['Label']=[1,2] 1245 labels = { 1246 'segmentation':hipp_bin, 1247 'description':hippleftORlabels, 1248 'HLProb':avgleft, 1249 'HRProb':avgright 1250 } 1251 return labels 1252 1253 1254def dap( x ): 1255 return( antspynet.deep_atropos( x, do_preprocessing=True, do_denoising=False )['segmentation_image'] ) 1256 1257def label_and_img_to_sr( img, label_img, sr_model, return_intensity=False, target_range=[1,0] ): 1258 """ 1259 Apply SR to a label image and the associated intensity. 1260 1261 img : the input intensity image 1262 1263 label_img : the input label image (segmentation) 1264 1265 sr_model : a super resolution model 1266 1267 return_intensity : boolean if True will return both intensity and labels, 1268 otherwise only the upsampled labels are returned. 1269 1270 target_range : eg [1,0] or [127.5,-127.5] are most likely 1271 1272 """ 1273 ulabs = np.unique( label_img.numpy() ) 1274 ulabs.sort() 1275 ulabs = ulabs[1:len(ulabs)] 1276 ulabs = ulabs.tolist() 1277 if return_intensity: 1278 return super_resolution_segmentation_per_label( 1279 img, 1280 label_img, 1281 [2,2,2], 1282 sr_model, 1283 segmentation_numbers=ulabs, 1284 target_range=target_range, 1285 max_lab_plus_one=True ) 1286 else: 1287 return super_resolution_segmentation_per_label( 1288 img, 1289 label_img, 1290 [2,2,2], 1291 sr_model, 1292 segmentation_numbers=ulabs, 1293 target_range=target_range, 1294 max_lab_plus_one=True )['super_resolution_segmentation'] 1295 1296def hierarchical_to_sr( t1hier, sr_model, tissue_sr=False, blending=0.5, verbose=False ): 1297 """ 1298 Apply SR to most output from the hierarchical function 1299 1300 t1hier : the hierarchical object 1301 1302 sr_model : a super resolution model 1303 1304 tissue_sr : boolean, if True will perform SR on the whole image which can 1305 be very memory intensive; only use if you have plenty of RAM. 1306 1307 blending : if not None, will blend the SR 1308 1309 verbose : boolean 1310 1311 """ 1312 # from >>> t1h.keys() 1313 img = t1hier['brain_n4_dnz'] 1314 myvarlist = [ 'mtl', 'cit168lab', 'snseg', 'bf', 'deep_cit168lab' ] 1315 for myvar in myvarlist: 1316 if verbose: 1317 print( myvar ) 1318 t1hier[myvar]=label_and_img_to_sr( img, t1hier[myvar], sr_model ) 1319 if verbose: 1320 print( 'dkt_cortex' ) 1321 temp = label_and_img_to_sr( img, t1hier['dkt_parc']['dkt_cortex'], sr_model, return_intensity=True ) 1322 tempupimg = temp['super_resolution'] 1323 if blending is not None: 1324 tempupimg = tempupimg * (1.0 - blending ) + ants.iMath( tempupimg, "Sharpen" ) * blending 1325 t1hier['dkt_parc']['dkt_cortex']= temp['super_resolution_segmentation'] 1326 if verbose: 1327 print( 'assemble summaries' ) 1328 t1hier['dataframes']["dktcortex"]= map_segmentation_to_dataframe( "dkt", t1hier['dkt_parc']['dkt_cortex'] ) 1329 t1hier['dataframes']["mtl"]=map_segmentation_to_dataframe( 'mtl_description', t1hier['mtl'] ) 1330 t1hier['dataframes']["cit168"]=map_segmentation_to_dataframe( 'CIT168_Reinf_Learn_v1_label_descriptions_pad', t1hier['cit168lab'] ) 1331 t1hier['dataframes']["snseg"]=map_segmentation_to_dataframe( 'CIT168_Reinf_Learn_v1_label_descriptions_pad', t1hier['snseg'] ) 1332 t1hier['dataframes']["bf"]=map_segmentation_to_dataframe( 'nbm3CH13', t1hier['bf'] ) 1333 t1hier['dataframes']["deep_cit168"]=map_segmentation_to_dataframe( 'CIT168_Reinf_Learn_v1_label_descriptions_pad', t1hier['deep_cit168lab'] ) 1334 1335 if tissue_sr: 1336 bmask = ants.threshold_image( t1hier['dkt_parc']['tissue_segmentation'], 1, 6 ) 1337 segcrop = ants.image_clone( t1hier['dkt_parc']['tissue_segmentation'] ) 1338 hemicrop = t1hier['left_right'] 1339 segcrop[ hemicrop == 2 ] = ( segcrop[ hemicrop == 2 ] + 6 ) 1340 segcrop = segcrop * bmask 1341 mysr = super_resolution_segmentation_per_label( 1342 t1hier['brain_n4_dnz'], segcrop, [2,2,2], sr_model, [1,2,3,4,5,6,7,8,9,10,11,12], 1343 dilation_amount=0, probability_images=None, 1344 probability_labels=[1,2,3,4,5,6,7,8,9,10,11,12], 1345 max_lab_plus_one=True, verbose=True ) 1346 if blending is not None: 1347 mysr['super_resolution'] = mysr['super_resolution'] * (1.0 - blending ) + ants.iMath( mysr['super_resolution'], "Sharpen" ) * blending 1348 # fix the doubled labels 1349 temp = mysr['super_resolution_segmentation'] 1350 for k in [7,8,9,10,11,12] : 1351 temp[ temp == k ] = temp[ temp == k ] - 6 1352 t1hier['brain_n4_dnz'] = mysr['super_resolution'] 1353 t1hier['dkt_parc']['tissue_segmentation'] = temp 1354 else: 1355 t1hier['brain_n4_dnz'] = tempupimg 1356 t1hier['dkt_parc']['tissue_segmentation'] = ants.resample_image_to_target( 1357 t1hier['dkt_parc']['tissue_segmentation'], tempupimg, 'nearestNeighbor') 1358 1359 if verbose: 1360 print("unique tissue labels") 1361 print( np.unique(t1hier['dkt_parc']['tissue_segmentation'].numpy()) ) 1362 1363 tissue = map_segmentation_to_dataframe( "tissues", 1364 ants.mask_image( 1365 t1hier['dkt_parc']['tissue_segmentation'], 1366 t1hier['dkt_parc']['tissue_segmentation'], [1,2,3,4,5,6] ) ) 1367 t1hier['dataframes']['tissues'] = tissue 1368 braintissuemask = ants.threshold_image( t1hier['dkt_parc']['tissue_segmentation'], 2, 6 ) 1369# 1370# mydcit = deep_cit168( t1hier['brain_n4_dnz'], binary_mask = braintissuemask ) 1371# t1hier['deep_cit168lab'] = mydcit['segmentation'] 1372# t1hier['dataframes']["deep_cit168"]=mydcit['description'] 1373 1374# deep_bf = deep_nbm( t1hier['brain_n4_dnz'] * braintissuemask, 1375# get_data("deep_nbm_rank",target_extension='.h5'), 1376# csfquantile=None, aged_template=True ) 1377# t1hier['bf'] = deep_bf['segmentation'] 1378# t1hier['dataframes']["bf"]=deep_bf['description'] 1379# 1380 return t1hier 1381 1382 1383def deep_mtl(t1, sr_model=None, verbose=True): 1384 1385 """ 1386 Hippocampal/Enthorhinal segmentation using "Deep Flash" 1387 1388 Perform hippocampal/entorhinal segmentation in T1 images using 1389 labels from Mike Yassa's lab 1390 1391 https://faculty.sites.uci.edu/myassa/ 1392 1393 The labeling is as follows: 1394 Label 0 : background 1395 Label 5 : left aLEC 1396 Label 6 : right aLEC 1397 Label 7 : left pMEC 1398 Label 8 : right pMEC 1399 Label 9 : left perirhinal 1400 Label 10: right perirhinal 1401 Label 11: left parahippocampal 1402 Label 12: right parahippocampal 1403 Label 13: left DG/CA3 1404 Label 14: right DG/CA3 1405 Label 15: left CA1 1406 Label 16: right CA1 1407 Label 17: left subiculum 1408 Label 18: right subiculum 1409 1410 """ 1411 1412 verbose = False 1413 1414 labels = (0, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18) 1415 label_descriptions = ['background', 1416 'left aLEC', 1417 'right aLEC', 1418 'left pMEC', 1419 'right pMEC', 1420 'left perirhinal', 1421 'right perirhinal', 1422 'left parahippocampal', 1423 'right parahippocampal', 1424 'left DG/CA3', 1425 'right DG/CA3', 1426 'left CA1', 1427 'right CA1', 1428 'left subiculum', 1429 'right subiculum' 1430 ] 1431 1432 template = ants.image_read(antspynet.get_antsxnet_data("deepFlashTemplateT1SkullStripped")) 1433 registration = ants.registration(fixed=template, moving=t1, 1434 type_of_transform="antsRegistrationSyNQuickRepro[a]", total_sigma=0.5, verbose=verbose) 1435 template_transforms = dict(fwdtransforms=registration['fwdtransforms'], 1436 invtransforms=registration['invtransforms']) 1437 t1_warped = registration['warpedmovout'] 1438 1439 df = antspynet.deep_flash(t1_warped, 1440 do_preprocessing=False, 1441 use_rank_intensity = True, 1442 verbose=verbose) 1443 1444 if False: 1445 if sr_model is not None: 1446 ulabs = np.unique( df['segmentation_image'].numpy() ) 1447 ulabs.sort() 1448 ulabs = ulabs[1:len(ulabs)] 1449 ulabs = ulabs.tolist() 1450 newprobs = super_resolution_segmentation_per_label( 1451 t1_warped, 1452 df['segmentation_image'], 1453 [2,2,2], 1454 sr_model, 1455 ulabs, 1456 probability_images=df['probability_images'][1:len(df['probability_images'])], 1457 probability_labels=ulabs, 1458 max_lab_plus_one=True ) 1459 for k in range(1, len(df['probability_images']) ): 1460 df['probability_images'][k] = newprobs['probability_images'][k-1] 1461 1462 if sr_model is not None: 1463 newprobs = super_resolution_segmentation_with_probabilities( 1464 t1_warped, df['probability_images'], sr_model ) 1465 df['probability_images'] = newprobs['sr_probabilities'] 1466 1467 if False: 1468 probability_images = list() 1469 for i in range(len(df['probability_images'])): 1470 probability_image = ants.apply_transforms(fixed=t1, 1471 moving=df['probability_images'][i], 1472 transformlist=template_transforms['invtransforms'], 1473 whichtoinvert=[True], 1474 interpolator="linear", 1475 verbose=verbose) 1476 probability_images.append(probability_image) 1477 image_matrix = ants.image_list_to_matrix(probability_images[1:(len(probability_images))], t1 * 0 + 1) 1478 background_foreground_matrix = np.stack([ants.image_list_to_matrix([probability_images[0]], t1 * 0 + 1), 1479 np.expand_dims(np.sum(image_matrix, axis=0), axis=0)]) 1480 foreground_matrix = np.argmax(background_foreground_matrix, axis=0) 1481 segmentation_matrix = (np.argmax(image_matrix, axis=0) + 1) * foreground_matrix 1482 segmentation_image = ants.matrix_to_images(np.expand_dims(segmentation_matrix, axis=0), t1 * 0 + 1)[0] 1483 relabeled_image = ants.image_clone(segmentation_image) 1484 for i in range(len(labels)): 1485 relabeled_image[segmentation_image==i] = labels[i] 1486 1487 relabeled_image = ants.apply_transforms( fixed=t1, 1488 moving=df['segmentation_image'], 1489 transformlist=template_transforms['invtransforms'], 1490 whichtoinvert=[True], 1491 interpolator="nearestNeighbor", 1492 verbose=verbose) 1493 mtl_description = map_segmentation_to_dataframe( 'mtl_description', relabeled_image ) 1494 1495 deep_mtl_dictionary = { 1496 'mtl_description':mtl_description, 1497 'mtl_segmentation':relabeled_image 1498# 'mtl_probability_images':probability_images 1499 } 1500 return(deep_mtl_dictionary) 1501 1502def localsyn( 1503 img, 1504 template, 1505 hemiS, 1506 templateHemi, 1507 whichHemi, 1508 padder, 1509 iterations, 1510 output_prefix, 1511 total_sigma=0.5 1512): 1513 """ 1514 Perform localized symmetric image registration between a subject hemisphere and a template. 1515 1516 This function extracts the specified hemisphere from both the subject and template images, 1517 identifies a tight region around the non-background portion of the template hemisphere, 1518 crops both images accordingly, and performs a symmetric normalization (SyN) registration 1519 between the cropped template and subject hemisphere. 1520 1521 Parameters 1522 ---------- 1523 img : ANTsImage 1524 The subject's anatomical image (e.g., T1-weighted) to be registered. 1525 1526 template : ANTsImage 1527 The template image to which the subject image will be registered. 1528 1529 hemiS : ANTsImage 1530 A hemisphere label image corresponding to the subject's space. Each hemisphere should 1531 be encoded with a distinct integer label (e.g., left = 1, right = 2). 1532 1533 templateHemi : ANTsImage 1534 A hemisphere label image in the template space, analogous to `hemiS`. 1535 1536 whichHemi : int 1537 The integer label for the hemisphere to process (e.g., 1 for left, 2 for right). 1538 1539 padder : int 1540 Number of morphological dilation steps (`MD`) applied to the cropped mask to ensure 1541 sufficient context around the hemisphere during cropping. 1542 1543 iterations : list of str 1544 Registration iteration schedule passed to `ants.registration()` (e.g., `["40x20x0"]`). 1545 1546 output_prefix : str 1547 Prefix to use for naming output files generated by the registration process. 1548 1549 total_sigma : float, optional 1550 The total Gaussian smoothing sigma used in the SyN registration (default: 0.5). 1551 1552 Returns 1553 ------- 1554 dict 1555 Dictionary containing output of `ants.registration()`, including the warped image, 1556 transforms, and potentially intermediate files. Keys include: 1557 - `'warpedmovout'`: warped subject image 1558 - `'fwdtransforms'`: forward transformation filenames 1559 - `'invtransforms'`: inverse transformation filenames 1560 - ... other ANTsPy registration outputs 1561 1562 Notes 1563 ----- 1564 - The function assumes that the subject and template hemisphere label images (`hemiS` and 1565 `templateHemi`) use the same labeling scheme and orientation. 1566 - The `padder` controls how much extra anatomical context is included around the cropped 1567 region. Larger values may improve registration robustness but increase computational cost. 1568 - The registration uses a reproducible variant of SyN (`antsRegistrationSyNQuickRepro[s]`) 1569 with a fixed random seed (1) for consistent outputs across runs. 1570 1571 Example 1572 ------- 1573 >>> syn_result = localsyn( 1574 ... img=subject_img, 1575 ... template=template_img, 1576 ... hemiS=subject_hemi_label, 1577 ... templateHemi=template_hemi_label, 1578 ... whichHemi=1, 1579 ... padder=10, 1580 ... iterations=["40x20x0"], 1581 ... output_prefix="subj_L_" 1582 ... ) 1583 >>> warped = syn_result['warpedmovout'] 1584 """ 1585 ihemi = img * ants.threshold_image(hemiS, whichHemi, whichHemi) 1586 themi = template * ants.threshold_image(templateHemi, whichHemi, whichHemi) 1587 loquant = np.quantile(themi.numpy(), 0.01) + 1e-6 # background threshold + epsilon 1588 hemicropmask = ants.threshold_image( 1589 templateHemi * ants.threshold_image(themi, loquant, math.inf), 1590 whichHemi, whichHemi 1591 ).iMath("MD", padder) 1592 tcrop = ants.crop_image(themi, hemicropmask) 1593 syn = ants.registration( 1594 tcrop, 1595 ihemi, 1596 'antsRegistrationSyNQuickRepro[s]', 1597 reg_iterations=iterations, 1598 flow_sigma=3.0, 1599 total_sigma=total_sigma, 1600 verbose=False, 1601 outprefix=output_prefix, 1602 random_seed=1 1603 ) 1604 return syn 1605 1606def hemi_reg( 1607 input_image, 1608 input_image_tissue_segmentation, 1609 input_image_hemisphere_segmentation, 1610 input_template, 1611 input_template_hemisphere_labels, 1612 output_prefix, 1613 padding=10, 1614 labels_to_register = [2,3,4,5], 1615 total_sigma=0.5, 1616 is_test=False ): 1617 """ 1618 hemisphere focused registration that will produce jacobians and figures to 1619 support data inspection 1620 1621 input_image: input image 1622 1623 input_image_tissue_segmentation: segmentation produced in ANTs style ie with 1624 labels defined by atropos brain segmentation (1 to 6) 1625 1626 input_image_hemisphere_segmentation: left (1) and right (2) hemisphere 1627 segmentation 1628 1629 input_template: template to which we register; prefer a population-specific 1630 relatively high-resolution template instead of MNI or biobank. 1631 1632 input_template_hemisphere_labels: a segmentation image of template hemispheres 1633 with left labeled 1 and right labeled 2 1634 1635 output_prefix: a path and prefix for registration related outputs 1636 1637 padding: number of voxels to pad images, needed for diffzero 1638 1639 labels_to_register: list of integer segmentation labels to use to define 1640 the tissue types / regions of the brain to register. 1641 1642 total_sigma: scalar >= 0.0 ; higher means more constrained registration. 1643 1644 is_test: boolean. this function can be long running by default. this would 1645 help testing more quickly by running fewer iterations. 1646 1647 """ 1648 1649 img = ants.rank_intensity( input_image ) 1650 # ionlycerebrum = ants.mask_image( input_image_tissue_segmentation, 1651 # input_image_tissue_segmentation, labels_to_register, 1 ) 1652 ionlycerebrum = brain_extraction( input_image ) 1653 1654 tonlycerebrum = brain_extraction( input_template ) 1655 template = ants.rank_intensity( input_template ) 1656 1657 regsegits=[200,200,20] 1658 1659# # upsample the template if we are passing SR as input 1660# if min(ants.get_spacing(img)) < 0.8 : 1661# regsegits=[200,200,200,20] 1662# template = ants.resample_image( template, (0.5,0.5,0.5), interp_type = 0 ) 1663# tonlycerebrum = ants.resample_image_to_target( tonlycerebrum, 1664# template, 1665# interp_type='nearestNeighbor', 1666# ) 1667 1668 if is_test: 1669 regsegits=[8,0,0] 1670 1671 input_template_hemisphere_labels = ants.resample_image_to_target( 1672 input_template_hemisphere_labels, 1673 template, 1674 interp_type='nearestNeighbor', 1675 ) 1676 1677 # now do a hemisphere focused registration 1678 synL = localsyn( 1679 img=img*ionlycerebrum, 1680 template=template*tonlycerebrum, 1681 hemiS=input_image_hemisphere_segmentation, 1682 templateHemi=input_template_hemisphere_labels, 1683 whichHemi=1, 1684 padder=padding, 1685 iterations=regsegits, 1686 output_prefix = output_prefix + "left_hemi_reg", 1687 total_sigma=total_sigma, 1688 ) 1689 synR = localsyn( 1690 img=img*ionlycerebrum, 1691 template=template*tonlycerebrum, 1692 hemiS=input_image_hemisphere_segmentation, 1693 templateHemi=input_template_hemisphere_labels, 1694 whichHemi=2, 1695 padder=padding, 1696 iterations=regsegits, 1697 output_prefix = output_prefix + "right_hemi_reg", 1698 total_sigma=total_sigma, 1699 ) 1700 1701 ants.image_write(synL['warpedmovout'], output_prefix + "left_hemi_reg.nii.gz" ) 1702 ants.image_write(synR['warpedmovout'], output_prefix + "right_hemi_reg.nii.gz" ) 1703 1704 fignameL = output_prefix + "_left_hemi_reg.png" 1705 ants.plot(synL['warpedmovout'],axis=2,ncol=8,nslices=24,filename=fignameL, black_bg=False, crop=True ) 1706 1707 fignameR = output_prefix + "_right_hemi_reg.png" 1708 ants.plot(synR['warpedmovout'],axis=2,ncol=8,nslices=24,filename=fignameR, black_bg=False, crop=True ) 1709 1710 lhjac = ants.create_jacobian_determinant_image( 1711 synL['warpedmovout'], 1712 synL['fwdtransforms'][0], 1713 do_log=1 1714 ) 1715 ants.image_write( lhjac, output_prefix+'left_hemi_jacobian.nii.gz' ) 1716 1717 rhjac = ants.create_jacobian_determinant_image( 1718 synR['warpedmovout'], 1719 synR['fwdtransforms'][0], 1720 do_log=1 1721 ) 1722 ants.image_write( rhjac, output_prefix+'right_hemi_jacobian.nii.gz' ) 1723 return { 1724 "synL":synL, 1725 "synLpng":fignameL, 1726 "synR":synR, 1727 "synRpng":fignameR, 1728 "lhjac":lhjac, 1729 "rhjac":rhjac 1730 } 1731 1732 1733 1734def region_reg( 1735 input_image, 1736 input_image_tissue_segmentation, 1737 input_image_region_segmentation, 1738 input_template, 1739 input_template_region_segmentation, 1740 output_prefix, 1741 padding=10, 1742 total_sigma=0.5, 1743 is_test=False ): 1744 """ 1745 region focused registration that will produce jacobians and figures to 1746 support data inspection. region-defining images should be binary. 1747 1748 input_image: input image 1749 1750 input_image_tissue_segmentation: segmentation produced in ANTs style ie with 1751 labels defined by atropos brain segmentation (1 to 6) 1752 1753 input_image_region_segmentation: a local region to register - binary. 1754 1755 input_template: template to which we register; prefer a population-specific 1756 relatively high-resolution template instead of MNI or biobank. brain extracted. 1757 1758 input_template_region_segmentation: a segmentation image of template regions - binary. 1759 1760 output_prefix: a path and prefix for registration related outputs 1761 1762 padding: number of voxels to pad images, needed for diffzero 1763 1764 total_sigma: scalar >= 0.0 ; higher means more constrained registration. 1765 1766 is_test: boolean. this function can be long running by default. this would 1767 help testing more quickly by running fewer iterations. 1768 1769 """ 1770 1771 img = ants.rank_intensity( input_image ) 1772 ionlycerebrum = brain_extraction( input_image ) 1773 template = ants.rank_intensity( input_template ) 1774 regsegits=[200,200,20] 1775 1776 # upsample the template if we are passing SR as input 1777 if min(ants.get_spacing(img)) < 0.8: 1778 regsegits=[200,200,200,20] 1779 template = ants.resample_image( template, (0.5,0.5,0.5), interp_type = 0 ) 1780 1781 if is_test: 1782 regsegits=[20,5,0] 1783 1784 input_template_region_segmentation = ants.resample_image_to_target( 1785 input_template_region_segmentation, 1786 template, 1787 interp_type='nearestNeighbor', 1788 ) 1789 1790 # now do a region focused registration 1791 synL = localsyn( 1792 img=img*ionlycerebrum, 1793 template=template, 1794 hemiS=input_image_region_segmentation, 1795 templateHemi=input_template_region_segmentation, 1796 whichHemi=1, 1797 padder=padding, 1798 iterations=regsegits, 1799 output_prefix = output_prefix + "region_reg", 1800 total_sigma=total_sigma, 1801 ) 1802 1803 ants.image_write(synL['warpedmovout'], output_prefix + "region_reg.nii.gz" ) 1804 1805 fignameL = output_prefix + "_region_reg.png" 1806 ants.plot(synL['warpedmovout'],axis=2,ncol=8,nslices=24,filename=fignameL, black_bg=False, crop=True ) 1807 1808 lhjac = ants.create_jacobian_determinant_image( 1809 synL['warpedmovout'], 1810 synL['fwdtransforms'][0], 1811 do_log=1 1812 ) 1813 ants.image_write( lhjac, output_prefix+'region_jacobian.nii.gz' ) 1814 1815 return { 1816 "synL":synL, 1817 "synLpng":fignameL, 1818 "lhjac":lhjac, 1819 'rankimg':img*ionlycerebrum, 1820 'template':template 1821 } 1822 1823 1824def t1_hypointensity( x, xsegmentation, xWMProbability, template, templateWMPrior, wmh_thresh=0.1 ): 1825 """ 1826 provide measurements that may help decide if a given t1 image is likely 1827 to have hypointensity. 1828 1829 x: input image; bias-corrected, brain-extracted and denoised 1830 1831 xsegmentation: input image hard-segmentation results 1832 1833 xWMProbability: input image WM probability 1834 1835 template: template image 1836 1837 templateWMPrior: template tissue prior 1838 1839 wmh_thresh: float used to threshold WMH probability and produce summary data 1840 1841 returns: 1842 1843 - wmh_summary: summary data frame based on thresholding WMH probability at wmh_thresh 1844 - wmh_probability_image: probability image denoting WMH probability; higher values indicate 1845 that WMH is more likely 1846 - wmh_evidence_of_existence: an integral evidence that indicates the likelihood that the input 1847 image content supports the presence of white matter hypointensity. 1848 greater than zero is supportive of WMH. the higher, the more so. 1849 less than zero is evidence against. 1850 - wmh_max_prob: max probability of wmh 1851 - features: the features driving WMH predictons 1852 1853 """ 1854 1855 if False: # Need to retrain this model with refined inference code 1856 return { 1857 "wmh_summary":None, 1858 "wmh_probability_image":None, 1859 "wmh_evidence_of_existence":None, 1860 "wmh_max_prob":None, 1861 "features":None } 1862 1863 mybig = [88,128,128] 1864 templatesmall = ants.resample_image( template, mybig, use_voxels=True ) 1865 qaff = ants.registration( 1866 ants.rank_intensity(x), 1867 ants.rank_intensity(templatesmall), 'antsRegistrationSyNQuickRepro[s]', 1868 syn_sampling=2, 1869 syn_metric='CC', 1870 reg_iterations = [25,15,0,0],total_sigma=0.5, 1871 aff_metric='GC', random_seed=1 ) 1872 afftx = qaff['fwdtransforms'][1] 1873 templateWMPrior2x = ants.apply_transforms( x, templateWMPrior, qaff['fwdtransforms'] ) 1874 cerebrum = ants.threshold_image( xsegmentation, 2, 4 ) 1875 realWM = ants.threshold_image( templateWMPrior2x , 0.1, math.inf ) 1876 inimg = ants.rank_intensity( x ) 1877 parcellateWMdnz = ants.kmeans_segmentation( inimg, 2, realWM, mrf=0.3 )['probabilityimages'][0] 1878 x2template = ants.apply_transforms( templatesmall, inimg, afftx, whichtoinvert=[True] ) 1879 parcellateWMdnz2template = ants.apply_transforms( templatesmall, 1880 cerebrum * parcellateWMdnz, afftx, whichtoinvert=[True] ) 1881 # features = rank+dnz-image, lprob, wprob, wprior at mybig resolution 1882 f1 = x2template.numpy() 1883 f2 = parcellateWMdnz2template.numpy() 1884 f3 = ants.apply_transforms( templatesmall, xWMProbability, afftx, whichtoinvert=[True] ).numpy() 1885 f4 = ants.apply_transforms( templatesmall, templateWMPrior, qaff['fwdtransforms'][0] ).numpy() 1886 myfeatures = np.stack( (f1,f2,f3,f4), axis=3 ) 1887 newshape = np.concatenate( [ [1],np.asarray( myfeatures.shape )] ) 1888 myfeatures = myfeatures.reshape( newshape ) 1889 1890 inshape = [None,None,None,4] 1891 wmhunet = antspynet.create_unet_model_3d( inshape, 1892 number_of_outputs = 1, 1893 number_of_layers = 4, 1894 mode = 'sigmoid' ) 1895 1896 wmhunet.load_weights( get_data("simwmhseg", target_extension='.h5') ) 1897 1898 pp = wmhunet.predict( myfeatures ) 1899 1900 limg = ants.from_numpy( tf.squeeze( pp[0] ).numpy( ) ) 1901 limg = ants.copy_image_info( templatesmall, limg ) 1902 lesresam = ants.apply_transforms( x, limg, afftx, whichtoinvert=[False] ) 1903 # lesresam = lesresam * cerebrum 1904 rnmdl = antspynet.create_resnet_model_3d( inshape, 1905 number_of_outputs = 1, 1906 layers = (1,2,3), 1907 residual_block_schedule = (3,4,6,3), squeeze_and_excite = True, 1908 lowest_resolution = 32, cardinality = 1, mode = "regression" ) 1909 rnmdl.load_weights( get_data("simwmdisc", target_extension='.h5' ) ) 1910 qq = rnmdl.predict( myfeatures ) 1911 1912 lesresamb = ants.threshold_image( lesresam, wmh_thresh, 1.0 ) 1913 lgo=ants.label_geometry_measures( lesresamb, lesresam ) 1914 wmhsummary = pd.read_csv( get_data("wmh_evidence", target_extension='.csv' ) ) 1915 wmhsummary.at[0,'Value']=lgo.at[0,'VolumeInMillimeters'] 1916 wmhsummary.at[1,'Value']=lgo.at[0,'IntegratedIntensity'] 1917 wmhsummary.at[2,'Value']=float(qq) 1918 1919 return { 1920 "wmh_summary":wmhsummary, 1921 "wmh_probability_image":lesresam, 1922 "wmh_evidence_of_existence":float(qq), 1923 "wmh_max_prob":lesresam.max(), 1924 "features":myfeatures } 1925 1926 1927 1928def deep_nbm( t1, 1929 nbm_weights, 1930 binary_mask=None, 1931 deform=False, 1932 aged_template=False, 1933 csfquantile=None, 1934 reflect=False, 1935 verbose=False ): 1936 1937 """ 1938 CH13 and Nucleus basalis of Meynert segmentation and subdivision 1939 1940 Perform CH13 and NBM segmentation in T1 images using Avants labels. 1941 1942 t1 : T1-weighted neuroimage antsImage - already brain extracted 1943 1944 nbm_weights : string weight file for parcellating unet 1945 1946 binary_mask : will restrict output to this mask 1947 1948 deform : boolean to correct for image deformation 1949 1950 aged_template : boolean to control which template to use 1951 1952 csfquantile: if not None, will try to remove CSF from the image. 1953 0.25 may be a good value to try. 1954 1955 verbose: boolean 1956 1957 The labeling is as follows: 1958 1959 Label,Description,Side 1960 1,CH13_left,left 1961 2,CH13_right,right 1962 3,NBM_left_ant,left 1963 4,NBM_left_mid,left 1964 5,NBM_left_pos,left 1965 6,NBM_right_ant,right 1966 7,NBM_right_mid,right 1967 8,NBM_right_pos,right 1968 1969 Failure modes will include odd image orientation (in which case you might 1970 use the registration option). A more nefarious issue can be a poor extraction 1971 of the cerebrum in the inferior frontal lobe. These can be unpredictable 1972 but if one sees a bad extraction, please check the mask that is output by 1973 this function to see if it excludes non-cerebral tissue. 1974 1975 """ 1976 1977 if not aged_template: 1978 refimg = ants.image_read( get_data( "CIT168_T1w_700um_pad", target_extension='.nii.gz' )) 1979 refimg = ants.rank_intensity( refimg ) 1980 refimg = ants.resample_image( refimg, [0.5,0.5,0.5] ) 1981 refimgseg = ants.image_read( get_data( "CIT168_basal_forebrain", target_extension='.nii.gz' )) 1982 refimgsmall = ants.resample_image( refimg, [2.0,2.0,2.0] ) 1983 else: 1984 refimg = ants.image_read( get_data( "CIT168_T1w_700um_pad_adni", target_extension='.nii.gz' )) 1985 refimg = ants.rank_intensity( refimg ) 1986 refimg = ants.resample_image( refimg, [0.5,0.5,0.5] ) 1987 refimgseg = ants.image_read( get_data( "CIT168_basal_forebrain_adni", target_extension='.nii.gz' )) 1988 refimgsmall = ants.resample_image( refimg, [2.0,2.0,2.0] ) 1989 1990 pt_labels = [1,2,3,4] 1991 group_labels = [0,1,2,3,4,5,6,7,8] 1992 reflection_labels = [0,2,1,6,7,8,3,4,5] 1993 crop_size = [144,96,64] 1994 1995 def nbmpreprocess( img, pt_labels, group_labels, masker=None, csfquantile=None, returndef=False, reflect=False ): 1996 1997 if masker is None: 1998 imgr = ants.rank_intensity( img ) 1999 else: 2000 imgr = ants.rank_intensity( img, mask = masker ) 2001 2002 if csfquantile is not None and masker is None: 2003 masker = ants.threshold_image( imgr, np.quantile(imgr[imgr>1e-4], csfquantile ), 1e9 ) 2004 2005 if masker is not None: 2006 imgr = imgr * masker 2007 2008 imgrsmall = ants.resample_image( imgr, [1,1,1] ) 2009 if not reflect: 2010 reg = ants.registration( 2011 refimgsmall, 2012 imgrsmall, 2013 'antsRegistrationSyNQuickRepro[s]', 2014 reg_iterations = [200,200,20],total_sigma=0.5, 2015 verbose=False ) 2016 else: 2017 myref = ants.reflect_image( imgrsmall, axis=0, tx='Translation' ) 2018 reg = ants.registration( 2019 refimgsmall, 2020 imgrsmall, 2021 'antsRegistrationSyNQuickRepro[s]', 2022 reg_iterations = [200,200,20], 2023 total_sigma=0.5, 2024 initial_transform = myref['fwdtransforms'][0], 2025 verbose=False ) 2026 if not returndef: 2027 imgraff = ants.apply_transforms( refimg, imgr, reg['fwdtransforms'][1], interpolator='linear' ) 2028 imgseg = ants.apply_transforms( refimg, refimgseg, reg['invtransforms'][1], interpolator='nearestNeighbor' ) 2029 else: 2030 imgraff = ants.apply_transforms( refimg, imgr, reg['fwdtransforms'], interpolator='linear' ) 2031 imgseg = ants.image_clone( refimgseg ) 2032 binseg = ants.mask_image( imgseg, imgseg, pt_labels, binarize=True ) 2033 imgseg = ants.mask_image( imgseg, imgseg, group_labels, binarize=False ) 2034 com = ants.get_center_of_mass( binseg ) 2035 return { 2036 "img": imgraff, 2037 "seg": imgseg, 2038 "imgc": special_crop( imgraff, com, crop_size ), 2039 "segc": special_crop( imgseg, com, crop_size ), 2040 "reg" : reg, 2041 "mask": masker 2042 } 2043 2044 2045 nLabels = len( group_labels ) 2046 number_of_outputs = len(group_labels) 2047 number_of_channels = 1 2048 ################################################ 2049 unet0 = antspynet.create_unet_model_3d( 2050 [ None, None, None, number_of_channels ], 2051 number_of_outputs = 1, # number of landmarks must be known 2052 number_of_layers = 4, # should optimize this wrt criterion 2053 number_of_filters_at_base_layer = 32, # should optimize this wrt criterion 2054 convolution_kernel_size = 3, # maybe should optimize this wrt criterion 2055 deconvolution_kernel_size = 2, 2056 pool_size = 2, 2057 strides = 2, 2058 dropout_rate = 0.0, 2059 weight_decay = 0, 2060 additional_options = "nnUnetActivationStyle", 2061 mode = "sigmoid" ) 2062 2063 unet1 = antspynet.create_unet_model_3d( 2064 [None,None,None,2], 2065 number_of_outputs=number_of_outputs, 2066 mode="classification", 2067 number_of_filters=(32, 64, 96, 128, 256), 2068 convolution_kernel_size=(3, 3, 3), 2069 deconvolution_kernel_size=(2, 2, 2), 2070 dropout_rate=0.0, 2071 weight_decay=0, 2072 additional_options = "nnUnetActivationStyle") 2073 2074 # concat output to input and pass to 2nd net 2075 # nextin = tf.concat( [ unet0.inputs[0], unet0.outputs[0] ], axis=4 ) 2076 import tensorflow as tf 2077 from tensorflow.keras.layers import Layer 2078 class myConcat(Layer): 2079 def call(self, x): 2080 return tf.concat(x, axis=4 ) 2081 nextin = myConcat()( [ unet0.inputs[0], unet0.outputs[0] ] ) 2082 2083 unetonnet = unet1( nextin ) 2084 unet_model = tf.keras.models.Model( 2085 unet0.inputs, 2086 [ unetonnet, unet0.outputs[0] ] ) 2087 2088 unet_model.load_weights( nbm_weights ) 2089 2090 imgprepro = nbmpreprocess( t1, pt_labels, group_labels, 2091 csfquantile = csfquantile, 2092 returndef = deform, 2093 reflect = reflect ) 2094 2095 def map_back( relo, t1, imgprepro, interpolator='linear', deform = False ): 2096 if not deform: 2097 relo = ants.apply_transforms( t1, relo, 2098 imgprepro['reg']['invtransforms'][0], whichtoinvert=[True], 2099 interpolator=interpolator ) 2100 else: 2101 relo = ants.apply_transforms( t1, relo, 2102 imgprepro['reg']['invtransforms'], interpolator=interpolator ) 2103 return relo 2104 2105 2106 #################################################### 2107 physspaceBF = imgprepro['imgc'] 2108 if verbose: 2109 print( "Mean preprocessed cropped image in nbm seg " + str(physspaceBF.mean() ) ) 2110 2111 tfarr1 = tf.cast( physspaceBF.numpy() ,'float32' ) 2112 newshapeBF = list( tfarr1.shape ) 2113 newshapeBF.insert(0,1) 2114 newshapeBF.insert(4,1) 2115 tfarr1 = tf.reshape(tfarr1, newshapeBF ) 2116 snpred = unet_model.predict( tfarr1 ) 2117 segpred = snpred[0] 2118 sigmoidpred = snpred[1] 2119 snpred1_image = ants.from_numpy( sigmoidpred[0,:,:,:,0] ) 2120 snpred1_image = ants.copy_image_info( physspaceBF, snpred1_image ) 2121 snpred1_image = map_back( snpred1_image, t1, imgprepro, 'linear', deform ) 2122 bint = ants.threshold_image( snpred1_image, 0.5, 1.0 ) 2123 probability_images = [] 2124 for jj in range(number_of_outputs-1): 2125 temp = ants.from_numpy( segpred[0,:,:,:,jj+1] ) 2126 temp = ants.copy_image_info( physspaceBF, temp ) 2127 temp = map_back( temp, t1, imgprepro, 'linear', deform ) 2128 probability_images.append( temp ) 2129 image_matrix = ants.image_list_to_matrix(probability_images, bint) 2130 segmentation_matrix = (np.argmax(image_matrix, axis=0) + 1) 2131 segmentation_image = ants.matrix_to_images(np.expand_dims(segmentation_matrix, axis=0), bint)[0] 2132 relabeled_image = ants.image_clone(segmentation_image) 2133 if not reflect: 2134 for i in range(1,len(group_labels)): 2135 relabeled_image[segmentation_image==(i)] = group_labels[i] 2136 else: 2137 for i in range(1,len(group_labels)): 2138 relabeled_image[segmentation_image==(i)] = reflection_labels[i] 2139 2140 bfsegdesc = map_segmentation_to_dataframe( 'nbm3CH13', relabeled_image ) 2141 2142 return { 'segmentation':relabeled_image, 2143 'description':bfsegdesc, 2144 'mask': imgprepro['mask'], 2145 'probability_images': probability_images } 2146 2147 2148def deep_nbm_old( t1, ch13_weights, nbm_weights, registration=True, 2149 csfquantile = 0.15, binary_mask=None, verbose=False ): 2150 2151 """ 2152 Nucleus basalis of Meynert segmentation and subdivision 2153 2154 Perform NBM segmentation in T1 images using Avants labels. 2155 2156 t1 : T1-weighted neuroimage antsImage - already brain extracted 2157 2158 ch13_weights : string weight file for ch13 2159 2160 nbm_weights : string weight file for nbm 2161 2162 registration : boolean to correct for image orientation and resolution by registration 2163 2164 csfquantile : float value below 0.5 that tries to trim residual CSF off brain. 2165 2166 binary_mask : will restrict output to this mask 2167 2168 The labeling is as follows: 2169 2170 Label,Description,Side 2171 1,CH13_left,left 2172 2,CH13_right,right 2173 3,NBM_left_ant,left 2174 4,NBM_left_mid,left 2175 5,NBM_left_pos,left 2176 6,NBM_right_ant,right 2177 7,NBM_right_mid,right 2178 8,NBM_right_pos,right 2179 2180 Failure modes will include odd image orientation (in which case you might 2181 use the registration option). A more nefarious issue can be a poor extraction 2182 of the cerebrum in the inferior frontal lobe. These can be unpredictable 2183 but if one sees a bad extraction, please check the mask that is output by 2184 this function to see if it excludes non-cerebral tissue. 2185 2186 """ 2187 2188 labels = [0, 1, 2, 3, 4, 5, 6, 7, 8] 2189 label_descriptions = ['background', 2190 'CH13_left', 2191 'CH13_right', 2192 'NBM_left_ant', 2193 'NBM_left_mid', 2194 'NBM_left_pos', 2195 'NBM_right_ant', 2196 'NBM_right_mid', 2197 'NBM_right_pos', 2198 ] 2199 2200 t1use = ants.iMath( t1, "Normalize" ) 2201 if registration: 2202 nbmtemplate = ants.image_read(get_data("nbm_template", target_extension=".nii.gz")) 2203 orireg = ants.registration( fixed = nbmtemplate, 2204 moving = t1use, 2205 type_of_transform="antsRegistrationSyNQuickRepro[a]", total_sigma=0.5,verbose=False ) 2206 t1use = orireg['warpedmovout'] 2207 2208 template = ants.image_read(get_data("CIT168_T1w_700um_pad_adni", target_extension=".nii.gz")) 2209 templateSmall = ants.resample_image( template, [2.0,2.0,2.0] ) 2210 registrationsyn = ants.registration( 2211 fixed=templateSmall, 2212 moving=ants.iMath(t1use,"Normalize"), 2213 type_of_transform="antsRegistrationSyNQuickRepro[s]", total_sigma=0.5, verbose=False ) 2214 2215 if verbose: 2216 print( registrationsyn['fwdtransforms'] ) 2217 2218 image = ants.iMath( t1use, "TruncateIntensity", 0.0001, 0.999 ).iMath("Normalize") 2219 bfPriorL1 = ants.image_read(get_data("CIT168_basal_forebrain_adni_prob_1_left", target_extension=".nii.gz")) 2220 bfPriorR1 = ants.image_read(get_data("CIT168_basal_forebrain_adni_prob_1_right", target_extension=".nii.gz")) 2221 bfPriorL2 = ants.image_read(get_data("CIT168_basal_forebrain_adni_prob_2_left", target_extension=".nii.gz")) 2222 bfPriorR2 = ants.image_read(get_data("CIT168_basal_forebrain_adni_prob_2_right", target_extension=".nii.gz")) 2223 2224 patchSize = [ 64, 64, 32 ] 2225 priorL1tosub = ants.apply_transforms( image, bfPriorL1, registrationsyn['invtransforms'] ).smooth_image( 3 ).iMath("Normalize") 2226 priorR1tosub = ants.apply_transforms( image, bfPriorR1, registrationsyn['invtransforms'] ).smooth_image( 3 ).iMath("Normalize") 2227 priorL2tosub = ants.apply_transforms( image, bfPriorL2, registrationsyn['invtransforms'] ).smooth_image( 3 ).iMath("Normalize") 2228 priorR2tosub = ants.apply_transforms( image, bfPriorR2, registrationsyn['invtransforms'] ).smooth_image( 3 ).iMath("Normalize") 2229 2230 if binary_mask is None: 2231 masker = ants.threshold_image(image, np.quantile(image[image>1e-4], csfquantile ), 1e9 ) 2232 else: 2233 masker = ants.apply_transforms( image, binary_mask, 2234 orireg['fwdtransforms'], interpolator='nearestNeighbor' ) 2235 2236 ch13point = ants.get_center_of_mass( priorL1tosub + priorR1tosub ) 2237 2238 nchanCH13 = 1 2239 nclasstosegCH13 = 3 # for ch13 2240 nchanNBM = 2 2241 nclasstosegNBM = 4 # for nbm 2242 actor = 'classification' 2243 nfilt = 32 2244 addoptsNBM = "nnUnetActivationStyle" 2245 unetCH13 = antspynet.create_unet_model_3d( 2246 [ None, None, None, nchanCH13 ], 2247 number_of_outputs = nclasstosegCH13, # number of landmarks must be known 2248 number_of_layers = 4, # should optimize this wrt criterion 2249 number_of_filters_at_base_layer = 32, # should optimize this wrt criterion 2250 convolution_kernel_size = 3, # maybe should optimize this wrt criterion 2251 deconvolution_kernel_size = 2, 2252 pool_size = 2, 2253 strides = 2, 2254 dropout_rate = 0, 2255 weight_decay = 0, 2256 mode = 'classification' ) 2257 unetCH13.load_weights( ch13_weights ) 2258 2259 physspace = special_crop( image, ch13point, patchSize) 2260 ch13array = physspace.numpy() 2261 newshape = list( ch13array.shape ) 2262 newshape.insert(0,1) 2263 newshape.append(1) 2264 ch13pred = unetCH13.predict( tf.reshape( ch13array, newshape ) ) 2265 probability_images = [] 2266 for jj in range(3): 2267 temp = ants.from_numpy( ch13pred[0,:,:,:,jj] ) 2268 probability_images.append( ants.copy_image_info( physspace, temp ) ) 2269 bint = physspace * 0 + 1 2270 image_matrix = ants.image_list_to_matrix(probability_images[1:(len(probability_images))], bint ) 2271 background_foreground_matrix = np.stack([ants.image_list_to_matrix([probability_images[0]], bint), 2272 np.expand_dims(np.sum(image_matrix, axis=0), axis=0)]) 2273 foreground_matrix = np.argmax(background_foreground_matrix, axis=0) 2274 segmentation_matrix = (np.argmax(image_matrix, axis=0) + 1) * foreground_matrix 2275 segmentation_image = ants.matrix_to_images(np.expand_dims(segmentation_matrix, axis=0), bint)[0] 2276 relabeled_image = ants.image_clone(segmentation_image) 2277 ch13totalback = ants.resample_image_to_target(relabeled_image, image, interp_type='nearestNeighbor') * masker 2278 if registration: 2279 ch13totalback = ants.apply_transforms( t1, ch13totalback, 2280 orireg['invtransforms'][0], whichtoinvert=[True], interpolator='nearestNeighbor' ) 2281 2282 if verbose: 2283 print("CH13 done") 2284 2285 maskind = 3 2286 nlayers = 4 # for unet 2287 unet1 = antspynet.create_unet_model_3d( 2288 [ None, None, None, 2 ], 2289 number_of_outputs = 1, # number of landmarks must be known 2290 number_of_layers = 4, # should optimize this wrt criterion 2291 number_of_filters_at_base_layer = 32, # should optimize this wrt criterion 2292 convolution_kernel_size = 3, # maybe should optimize this wrt criterion 2293 deconvolution_kernel_size = 2, 2294 pool_size = 2, 2295 strides = 2, 2296 dropout_rate = 0.0, 2297 weight_decay = 0, 2298 additional_options = "nnUnetActivationStyle", 2299 mode = "sigmoid" ) 2300 maskinput = tf.keras.layers.Input( [ None, None, None, 1 ] ) 2301 posteriorMask1 = tf.keras.layers.multiply( 2302 [ unet1.outputs[0] , maskinput ], name='maskTimesPosteriors1' ) 2303 unet = tf.keras.models.Model( [ unet1.inputs[0], maskinput ], posteriorMask1 ) 2304 2305 unet2 = antspynet.create_unet_model_3d( 2306 [ None, None, None, 2 ], 2307 number_of_outputs = nclasstosegNBM, # number of landmarks must be known 2308 number_of_layers = 4, # should optimize this wrt criterion 2309 number_of_filters_at_base_layer = 32, # should optimize this wrt criterion 2310 convolution_kernel_size = 3, # maybe should optimize this wrt criterion 2311 deconvolution_kernel_size = 2, 2312 pool_size = 2, 2313 strides = 2, 2314 dropout_rate = 0.0, 2315 weight_decay = 0, 2316 additional_options = "nnUnetActivationStyle", 2317 mode = "classification" ) 2318 2319 temp = tf.split( unet1.inputs[0], 2, axis=4 ) 2320 temp[1] = unet.outputs[0] 2321 newmult = tf.concat( temp, axis=4 ) 2322 unetonnet = unet2( newmult ) 2323 unetNBM = tf.keras.models.Model( 2324 unet.inputs, 2325 [ unetonnet, unet.outputs[0] ] ) 2326 unetNBM.load_weights( nbm_weights ) 2327 2328 # do each side separately 2329 bfseg = t1 * 0.0 2330 for nbmnum in [0,1]: 2331 if nbmnum == 0: 2332 nbmpoint = ants.get_center_of_mass( priorL2tosub ) 2333 nbmprior = special_crop( priorL2tosub, nbmpoint, patchSize).numpy() # prior 2334 labels=[3,4,5] 2335 if nbmnum == 1: 2336 nbmpoint = ants.get_center_of_mass( priorR2tosub ) 2337 nbmprior = special_crop( priorR2tosub, nbmpoint, patchSize).numpy() # prior 2338 labels=[6,7,8] 2339 physspaceNBM = special_crop( image, nbmpoint, patchSize) # image 2340 nbmmask = special_crop( masker, nbmpoint, patchSize).numpy() # mask 2341 tfarr1 = tf.stack( [physspaceNBM.numpy(),nbmprior], axis=3 ) 2342 newshapeNBM = list( tfarr1.shape ) 2343 newshapeNBM.insert(0,1) 2344 tfarr1 = tf.reshape(tfarr1, newshapeNBM ) 2345 tfarr2 = tf.reshape( nbmmask, newshape ) 2346 nbmpred = unetNBM.predict( ( tfarr1, tfarr2 ) ) 2347 segpred = nbmpred[0] 2348 sigmoidpred = nbmpred[1] 2349 nbmpred1_image = ants.from_numpy( sigmoidpred[0,:,:,:,0] ) 2350 nbmpred1_image = ants.copy_image_info( physspaceNBM, nbmpred1_image ) 2351 bint = ants.threshold_image( nbmpred1_image, 0.5, 1.0 ).iMath("GetLargestComponent",1) 2352 probability_images = [] 2353 for jj in range(3): 2354 temp = ants.from_numpy( segpred[0,:,:,:,jj+1] ) 2355 probability_images.append( ants.copy_image_info( physspaceNBM, temp ) ) 2356 image_matrix = ants.image_list_to_matrix(probability_images, bint) 2357 segmentation_matrix = (np.argmax(image_matrix, axis=0) + 1) 2358 segmentation_image = ants.matrix_to_images(np.expand_dims(segmentation_matrix, axis=0), bint)[0] 2359 relabeled_image = ants.image_clone(segmentation_image) 2360 for i in range(len(labels)): 2361 relabeled_image[segmentation_image==(i+1)] = labels[i] 2362 relabeled_image = ants.resample_image_to_target(relabeled_image, image, interp_type='nearestNeighbor') 2363 if registration: 2364 relabeled_image = ants.apply_transforms( t1, relabeled_image, 2365 orireg['invtransforms'][0], whichtoinvert=[True], 2366 interpolator='nearestNeighbor' ) 2367 if verbose: 2368 print("NBM" + str( nbmnum ) ) 2369 bfseg = bfseg + relabeled_image 2370 bfseg = ch13totalback + bfseg * ants.threshold_image( ch13totalback, 0, 0 ) 2371 bfsegdesc = map_segmentation_to_dataframe( 'nbm3CH13', bfseg ) 2372 2373 if registration: 2374 masker = ants.apply_transforms( t1, masker, 2375 orireg['invtransforms'][0], whichtoinvert=[True], 2376 interpolator='nearestNeighbor' ) 2377 2378 return { 'segmentation':bfseg, 'description':bfsegdesc, 'mask': masker } 2379 2380 2381 2382 2383 2384def deep_cit168( t1, binary_mask = None, 2385 syn_type='antsRegistrationSyNQuickRepro[s]', 2386 priors = None, verbose = False): 2387 2388 """ 2389 CIT168 atlas segmentation with a parcellation unet. 2390 2391 Perform CIT168 segmentation in T1 images using Pauli atlas (CIT168) labels. 2392 2393 t1 : T1-weighted neuroimage antsImage - already brain extracted. image should 2394 be normalized 0 to 1 and with a nicely uniform histogram (no major outliers). 2395 we do a little work internally to help this but no guarantees it will handle 2396 all possible confounding issues. 2397 2398 binary_mask : will restrict output to this mask 2399 2400 syn_type : the type of registration used for generating priors; usually 2401 either SyN or antsRegistrationSyNQuickRepro[s] for repeatable results 2402 2403 priors : the user can provide their own priors through this argument; for 2404 example, the user may run this function twice, with the output of the first 2405 giving input to the second run. 2406 2407 verbose: boolean 2408 2409 Failure modes will primarily occur around red nucleus and caudate nucleus. 2410 For the latter, one might consider masking by the ventricular CSF, in particular 2411 near the anterior and inferior portion of the caudate in subjects with 2412 large ventricles. Low quality images with high atropy are also likely outside 2413 of the current range of the trained models. Iterating the model may help. 2414 2415 """ 2416 def tfsubset( x, indices ): 2417 with tf.device('/CPU:0'): 2418 outlist = [] 2419 for k in indices: 2420 outlist.append( x[:,:,:,int(k)] ) 2421 return tf.stack( outlist, axis=3 ) 2422 2423 def tfsubsetbatch( x, indices ): 2424 with tf.device('/CPU:0'): 2425 outlist2 = [] 2426 for j in range( len( x ) ): 2427 outlist = [] 2428 for k in indices: 2429 if len( x[j].shape ) == 5: 2430 outlist.append( x[j][k,:,:,:,:] ) 2431 if len( x[j].shape ) == 4: 2432 outlist.append( x[j][k,:,:,:] ) 2433 outlist2.append( tf.stack( outlist, axis=0 ) ) 2434 return outlist2 2435 2436 2437 registration = True 2438 cit168seg = t1 * 0 2439 myprior = ants.image_read(get_data("det_atlas_25_pad_LR_adni", target_extension=".nii.gz")) 2440 nbmtemplate = ants.image_read( get_data( "CIT168_T1w_700um_pad_adni", target_extension=".nii.gz" ) ) 2441 nbmtemplate = ants.resample_image( nbmtemplate, [0.5,0.5,0.5] ) 2442 templateSmall = ants.resample_image( nbmtemplate, [2.0,2.0,2.0] ) 2443 orireg = ants.registration( 2444 fixed = templateSmall, 2445 moving = ants.iMath( t1, "Normalize" ), 2446 type_of_transform=syn_type, total_sigma=0.5, verbose=False ) 2447 image = ants.apply_transforms( nbmtemplate, ants.iMath( t1, "Normalize" ), 2448 orireg['fwdtransforms'][1] ) 2449 image = ants.iMath( image, "TruncateIntensity",0.001,0.999).iMath("Normalize") 2450 patchSize = [ 160,160,112 ] 2451 if priors is None: 2452 priortosub = ants.apply_transforms( image, myprior, 2453 orireg['invtransforms'][1], interpolator='nearestNeighbor' ) 2454 else: 2455 if verbose: 2456 print("using priors") 2457 priortosub = ants.apply_transforms( image, priors, 2458 orireg['fwdtransforms'][1], interpolator='nearestNeighbor' ) 2459 bmask = ants.threshold_image( priortosub, 1, 999 ) 2460 # this identifies the cropping location - assumes a good registration 2461 pt = list( ants.get_center_of_mass( bmask ) ) 2462 pt[1] = pt[1] + 10.0 # same as we did in training 2463 2464 physspaceCIT = special_crop( image, pt, patchSize) # image 2465 2466 if binary_mask is not None: 2467 binary_mask_use = ants.apply_transforms( nbmtemplate, binary_mask, 2468 orireg['fwdtransforms'][1] ) 2469 binary_mask_use = special_crop( binary_mask_use, pt, patchSize) 2470 2471 for sn in [True,False]: 2472 if sn: 2473 group_labels = [0,7,8,9,23,24,25,33,34] 2474 newfn=get_data( "deepCIT168_sn", target_extension=".h5" ) 2475 else: 2476 group_labels = [0,1,2,5,6,17,18,21,22] 2477 newfn=get_data( "deepCIT168", target_extension=".h5" ) 2478 2479 number_of_outputs = len(group_labels) 2480 number_of_channels = len(group_labels) 2481 2482 unet0 = antspynet.create_unet_model_3d( 2483 [ None, None, None, number_of_channels ], 2484 number_of_outputs = 1, # number of landmarks must be known 2485 number_of_layers = 4, # should optimize this wrt criterion 2486 number_of_filters_at_base_layer = 32, # should optimize this wrt criterion 2487 convolution_kernel_size = 3, # maybe should optimize this wrt criterion 2488 deconvolution_kernel_size = 2, 2489 pool_size = 2, 2490 strides = 2, 2491 dropout_rate = 0.0, 2492 weight_decay = 0, 2493 additional_options = "nnUnetActivationStyle", 2494 mode = "sigmoid" ) 2495 2496 unet1 = antspynet.create_unet_model_3d( 2497 [None,None,None,2], 2498 number_of_outputs=number_of_outputs, 2499 mode="classification", 2500 number_of_filters=(32, 64, 96, 128, 256), 2501 convolution_kernel_size=(3, 3, 3), 2502 deconvolution_kernel_size=(2, 2, 2), 2503 dropout_rate=0.0, 2504 weight_decay=0, 2505 additional_options = "nnUnetActivationStyle") 2506 2507# temp = tf.split( unet0.inputs[0], 9, axis=4 ) 2508 import tensorflow as tf 2509 from tensorflow.keras.layers import Layer 2510 class mySplit(Layer): 2511 def call(self, x): 2512 return tf.split(x, 9, axis=4 ) 2513 temp = mySplit()( unet0.inputs[0] ) 2514 2515 temp[1] = unet0.outputs[0] 2516 2517 class myConcat(Layer): 2518 def call(self, x): 2519 return tf.concat(x, axis=4 ) 2520# newmult = tf.concat( temp[0:2], axis=4 ) 2521 newmult = myConcat()( temp[0:2] ) 2522 unetonnet = unet1( newmult ) 2523 unet_model = tf.keras.models.Model( 2524 unet0.inputs, 2525 [ unetonnet, unet0.outputs[0] ] ) 2526 unet_model.load_weights( newfn ) 2527 ################### 2528 nbmprior = special_crop( priortosub, pt, patchSize).numpy() # prior 2529 imgnp = tf.reshape( physspaceCIT.numpy(), [160, 160, 112,1]) 2530 nbmprior = tf.one_hot( nbmprior, 35 ) 2531 nbmprior = tfsubset( nbmprior, group_labels[1:len(group_labels)] ) 2532 imgnp = tf.reshape( tf.concat( [imgnp,nbmprior], axis=3), [1,160, 160, 112,9]) 2533 2534 nbmpred = unet_model.predict( imgnp ) 2535 segpred = nbmpred[0] 2536 sigmoidpred = nbmpred[1] 2537 2538 def map_back_cit( relo, t1, orireg, interpolator='linear' ): 2539 relo = ants.apply_transforms( t1, relo, 2540 orireg['invtransforms'][0], whichtoinvert=[True], 2541 interpolator=interpolator ) 2542 return relo 2543 2544 nbmpred1_image = ants.from_numpy( sigmoidpred[0,:,:,:,0] ) 2545 nbmpred1_image = ants.copy_image_info( physspaceCIT, nbmpred1_image ) 2546 nbmpred1_image = map_back_cit( nbmpred1_image, t1, orireg, 'linear' ) 2547 if binary_mask is not None: 2548 nbmpred1_image = nbmpred1_image * binary_mask 2549 bint = ants.threshold_image( nbmpred1_image, 0.5, 1.0 ) 2550 2551 probability_images = [] 2552 for jj in range(1,len(group_labels)): 2553 temp = ants.from_numpy( segpred[0,:,:,:,jj] ) 2554 temp = ants.copy_image_info( physspaceCIT, temp ) 2555 temp = map_back_cit( temp, t1, orireg, 'linear' ) 2556 probability_images.append( temp ) 2557 2558 image_matrix = ants.image_list_to_matrix(probability_images, bint) 2559 segmentation_matrix = (np.argmax(image_matrix, axis=0) + 1) 2560 segmentation_image = ants.matrix_to_images(np.expand_dims(segmentation_matrix, axis=0), bint)[0] 2561 relabeled_image = ants.image_clone(segmentation_image*0.) 2562 for i in np.unique(segmentation_image.numpy()): 2563 if i > 0 : 2564 temp = ants.threshold_image(segmentation_image,i,i) 2565 if group_labels[int(i)] < 33: 2566 temp = ants.iMath( temp, "GetLargestComponent",1) 2567 relabeled_image = relabeled_image + temp*group_labels[int(i)] 2568 cit168seg = cit168seg + relabeled_image 2569 2570 cit168segdesc = map_segmentation_to_dataframe( 'CIT168_Reinf_Learn_v1_label_descriptions_pad', cit168seg ).dropna(axis=0) 2571 2572 return { 'segmentation':cit168seg, 'description':cit168segdesc } 2573 2574 2575def preprocess_intensity( x, brain_extraction, 2576 intensity_truncation_quantiles=[1e-4, 0.999], 2577 rescale_intensities=True ): 2578 """ 2579 Default intensity processing for a brain-extracted T1-weighted image. 2580 2581 Arguments 2582 --------- 2583 2584 x : T1-weighted neuroimage antsImage after brain extraction applied 2585 2586 brain_extraction : T1-weighted neuroimage brain extraction / segmentation 2587 2588 intensity_truncation_quantiles: parameters passed to TruncateIntensity; the 2589 first value truncates values below this quantile; the second truncates 2590 values above this quantile. 2591 2592 rescale_intensities: boolean passed to n4 2593 2594 Returns 2595 ------- 2596 processed image 2597 """ 2598 brain_extraction = ants.resample_image_to_target( brain_extraction, x, interp_type='nearestNeighbor' ) 2599 img = x * brain_extraction 2600 # print( "PP input : " + str( img.mean() ) ) 2601 img = ants.iMath( img, "TruncateIntensity", intensity_truncation_quantiles[0], intensity_truncation_quantiles[1] ).iMath( "Normalize" ) 2602 # print( "PP Trunc : " + str( img.mean() ) ) 2603 # img = ants.denoise_image( img, brain_extraction, noise_model='Gaussian') 2604 # print( "PP DNZ : " + str( img.mean() ) ) 2605 img = ants.n4_bias_field_correction( img, mask=brain_extraction, rescale_intensities=rescale_intensities, ).iMath("Normalize") 2606 # print( "PP N4 : " + str( img.mean() ) ) 2607 return img 2608 2609 2610def hierarchical( x, output_prefix, labels_to_register=[2,3,4,5], 2611 imgbxt=None, img6seg=None, cit168 = False, is_test=False, 2612 atropos_prior=None, 2613 sr_model = None, 2614 verbose=True ): 2615 """ 2616 Default processing for a T1-weighted image. See README. 2617 2618 Arguments 2619 --------- 2620 x : T1-weighted neuroimage antsImage 2621 2622 output_prefix: string directory and prefix 2623 2624 labels_to_register: list of integer segmentation labels (of 1 to 6 as defined 2625 by atropos: csf, gm, wm, dgm, brainstem, cerebellum) to define 2626 the tissue types / regions of the brain to register. set to None to 2627 skip registration which will be faster but omit some results. 2628 2629 imgbxt : pre-existing brain extraction - a binary image - will disable some processing 2630 2631 img6seg: pre-existing brain segmentation - a 6-class image ( ANTs standard ) 2632 2633 cit168 : boolean returns labels from CIT168 atlas with high-resolution registration 2634 otherwise, low-resolution regitration is used. 2635 2636 is_test: boolean ( parameters to run more quickly but with low quality ) 2637 2638 atropos_prior: prior weight for atropos post-processing; set to None if you 2639 do not want to use this. will modify the CSF, GM and WM segmentation to 2640 better fit the image intensity at the resolution of the input image. 2641 2642 sr_model: optional will trigger sr-based upsampling in some functions. 2643 2644 verbose: boolean 2645 2646 Returns 2647 ------- 2648 dataframes and associated derived data 2649 2650 - brain_n4_dnz : extracted brain denoised and bias corrected 2651 - brain_extraction : brain mask 2652 - rbp: random basis projection results 2653 - left_right : left righ hemisphere segmentation 2654 - dkt_parc : dictionary object containing segmentation labels 2655 - registration : dictionary object containing registration results 2656 - hippLR : dictionary object containing hippocampus results 2657 - medial_temporal_lobe : dictionary object containing deep_flash (medial temporal lobe parcellation) results 2658 - white_matter_hypointensity : dictionary object containing WMH results 2659 - wm_tractsL : white matter tracts, left 2660 - wm_tractsR : white matter tracts, right 2661 - dataframes : summary data frames 2662 2663 """ 2664 if x.dimension != 3: 2665 raise ValueError('hierarchical: input image should be 3-dimensional') 2666 2667 if verbose: 2668 print("Read") 2669 tfn = get_data('T_template0', target_extension='.nii.gz' ) 2670 tfnw = get_data('T_template0_WMP', target_extension='.nii.gz' ) 2671 tlrfn = get_data('T_template0_LR', target_extension='.nii.gz' ) 2672 bfn = antspynet.get_antsxnet_data( "croppedMni152" ) 2673 2674 ##### read images and do simple bxt ops 2675 templatea = ants.image_read( tfn ) 2676 if verbose: 2677 print("bxt") 2678 templatea = ( templatea * antspynet.brain_extraction( templatea, 't1' ) ).iMath( "Normalize" ) 2679 templateawmprior = ants.image_read( tfnw ) 2680 templatealr = ants.image_read( tlrfn ) 2681 templateb = ants.image_read( bfn ) 2682 templateb = ( templateb * antspynet.brain_extraction( templateb, 't1' ) ).iMath( "Normalize" ) 2683 if imgbxt is None: 2684 probablySR = False 2685 # brain_extraction( ants.iMath( x, "Normalize" ) ) 2686 imgbxt = antspynet.brain_extraction( ants.iMath( x, "Normalize" ), modality="t1threetissue")['segmentation_image'].threshold_image(1,1) 2687 img = preprocess_intensity( ants.iMath( x, "Normalize" ), imgbxt ) 2688 else: 2689 probablySR = True 2690 img = ants.iMath( x, "Normalize" ) 2691 2692 if verbose: 2693 print("QC -- mean of img : "+str(img.mean())) 2694 print(". -- mean of bxt : "+str(imgbxt.mean())) 2695 2696 # this is an unbiased method for identifying predictors that can be used to 2697 # rank / sort data into clusters, some of which may be associated 2698 # with outlierness or low-quality data 2699 templatesmall = ants.resample_image( templateb, (91,109,91), use_voxels=True ) 2700 myqc = inspect_raw_t1( img, output_prefix=output_prefix, option='brain' ) 2701 2702 if verbose: 2703 print("qc results") 2704 print( myqc ) 2705 2706 ##### intensity modifications 2707 img = ants.iMath( img, "Normalize" ) 2708 2709 if verbose: 2710 print("parcellation -- mean of img " + str(img.mean())) 2711 2712 ##### hierarchical labeling 2713 myparc = deep_brain_parcellation( x, templateb, 2714 img6seg = img6seg, 2715 atropos_prior = atropos_prior, 2716 do_cortical_propagation = not is_test, verbose=verbose ) 2717 2718 2719 cit168lab = None 2720 cit168reg = None 2721 cit168lab_desc = None 2722 cit168adni = get_data( "CIT168_T1w_700um_pad_adni",target_extension='.nii.gz') 2723 cit168adni = ants.image_read( cit168adni ).iMath("Normalize") 2724 cit168labT = get_data( "det_atlas_25_pad_LR_adni", target_extension='.nii.gz' ) 2725 cit168labT = ants.image_read( cit168labT ) 2726 cit168labStem = get_data( "CIT168_T1w_700um_pad_adni_brainstem", target_extension='.nii.gz' ) 2727 cit168labStem = ants.image_read( cit168labStem ) 2728 2729 if verbose: 2730 print("cit168-- mean of input image : " + str(img.mean())) 2731 print("cit168-- mean of input imgbxt : " + str(imgbxt.mean())) 2732 print(". -- mean of input segmentation : " + str(myparc['tissue_segmentation'].mean())) 2733 2734 cit168reg = region_reg( 2735 input_image = img, 2736 input_image_tissue_segmentation=myparc['tissue_segmentation'], 2737 input_image_region_segmentation=imgbxt, 2738 input_template=cit168adni, 2739 input_template_region_segmentation=ants.threshold_image( cit168adni, 0.15, 1 ), 2740 output_prefix=output_prefix + "_CIT168RRSYN", 2741 padding=10, 2742 total_sigma=0.5, 2743 is_test=not cit168 )['synL'] 2744 cit168lab = ants.apply_transforms( img, cit168labT, 2745 cit168reg['invtransforms'], interpolator = 'nearestNeighbor' ) 2746 cit168lab_desc = map_segmentation_to_dataframe( 'CIT168_Reinf_Learn_v1_label_descriptions_pad', cit168lab ).dropna(axis=0) 2747 2748 # optional - quick look at result 2749 # bxt_png = output_prefix + "_brain_extraction_dnz_n4_view.png" 2750 # ants.plot(img * 255.0,axis=2,ncol=8,nslices=24, crop=True, black_bg=False, 2751 # filename = bxt_png ) 2752 2753 if verbose: 2754 print("hemi -- mean of input image : " + str(img.mean())) 2755 2756 # assuming data is reasonable quality, we should proceed with the rest ... 2757 mylr = label_hemispheres( img, templatea, templatealr ) 2758 2759 ##### accumulate data into data frames 2760 hemi = map_segmentation_to_dataframe( "hemisphere", myparc['hemisphere_labels'] ) 2761 tissue = map_segmentation_to_dataframe( "tissues", myparc['tissue_segmentation'] ) 2762 dktl = map_segmentation_to_dataframe( "lobes", myparc['dkt_lobes'] ) 2763 dktp = map_segmentation_to_dataframe( "dkt", myparc['dkt_parcellation'] ) 2764 dktc = None 2765 if not is_test: 2766 dktc = map_segmentation_to_dataframe( "dkt", myparc['dkt_cortex'] ) 2767 2768 tissue_seg_png = output_prefix + "_seg.png" 2769 ants.plot( img*255, myparc['tissue_segmentation'], axis=2, nslices=21, ncol=7, 2770 alpha=0.6, filename=tissue_seg_png, 2771 crop=True, black_bg=False ) 2772 2773 myhypo = None 2774# if verbose: 2775# print("WMH") 2776 ##### below here are more exploratory nice to have outputs 2777# myhypo = t1_hypointensity( 2778# img, 2779# myparc['tissue_segmentation'], # segmentation 2780# myparc['tissue_probabilities'][3], # wm posteriors 2781# templatea, 2782# templateawmprior ) 2783 2784 if verbose: 2785 print("icv") 2786 2787 myicv = icv( x ) 2788 2789 braintissuemask = ants.threshold_image( myparc['tissue_segmentation'], 2, 6 ) 2790 if verbose: 2791 print( myicv ) 2792 print("NBM : mean of input image " + str( (img * braintissuemask).mean() ) ) 2793 ##### deep_nbm basal forebrain parcellation 2794 deep_bf = deep_nbm( img * braintissuemask, 2795 get_data("deep_nbm_rank",target_extension='.h5'), 2796 csfquantile=None, aged_template=True, verbose=verbose ) 2797 2798 if is_test: 2799 mydataframes = { 2800 "icv": myicv, 2801 "rbp": myqc['brain'], 2802 "tissues":tissue, 2803 "dktlobes":dktl, 2804 "dktregions":dktp, 2805 "dktcortex":dktc, 2806 "bf":deep_bf['description'], 2807 } 2808 2809 outputs = { 2810 "brain_n4_dnz": img, 2811 "brain_extraction": imgbxt, 2812 "left_right": mylr, 2813 "dkt_parc": myparc, 2814 "bf":deep_bf['segmentation'], 2815 "dataframes": mydataframes 2816 } 2817 2818 return outputs 2819 2820 2821 ##### traditional deformable registration as a high-resolution complement to above 2822 wm_tractsL = None 2823 wm_tractsR = None 2824 wmtdfL = None 2825 wmtdfR = None 2826 reg = None 2827 if labels_to_register is not None and not is_test: 2828 reg = hemi_reg( 2829 input_image = img, 2830 input_image_tissue_segmentation = myparc['tissue_segmentation'], 2831 input_image_hemisphere_segmentation = mylr, 2832 input_template=templatea, 2833 input_template_hemisphere_labels=templatealr, 2834 output_prefix = output_prefix + "_SYN", 2835 labels_to_register = labels_to_register, 2836 is_test=is_test ) 2837 if verbose: 2838 print("wm tracts") 2839 ##### how to use the hemi-reg output to generate any roi value from a template roi 2840 wm_tracts = ants.image_read( get_data( "wm_major_tracts", target_extension='.nii.gz' ) ) 2841 wm_tractsL = ants.apply_transforms( img, wm_tracts, reg['synL']['invtransforms'], 2842 interpolator='nearestNeighbor' ) * ants.threshold_image( mylr, 1, 1 ) 2843 wm_tractsR = ants.apply_transforms( img, wm_tracts, reg['synR']['invtransforms'], 2844 interpolator='nearestNeighbor' ) * ants.threshold_image( mylr, 2, 2 ) 2845 wmtdfL = map_segmentation_to_dataframe( "wm_major_tracts", wm_tractsL ) 2846 wmtdfR = map_segmentation_to_dataframe( "wm_major_tracts", wm_tractsR ) 2847 2848 if verbose: 2849 print("hippocampus") 2850 2851 ##### specialized labeling for hippocampus 2852 do_hippmapper = False 2853 if do_hippmapper: 2854 ntries = 10 2855 if is_test: 2856 ntries = 1 2857 hippLR = deep_hippo( img=img, template=templateb, number_of_tries=ntries, 2858 tx_type='Affine', sr_model = sr_model ) 2859 2860 if verbose: 2861 print("medial temporal lobe") 2862 2863 ##### deep_flash medial temporal lobe parcellation 2864 deep_flash = deep_mtl(img, sr_model = sr_model ) 2865 2866 if verbose: 2867 print("deep CIT168") 2868 ##### deep CIT168 segmentation - relatively fast 2869 deep_cit = deep_cit168( img, priors = cit168lab, 2870 binary_mask = braintissuemask ) 2871 2872 if verbose: 2873 print( "SN-specific segmentation" ) 2874# input_image_region_segmentation, input_template, input_template_region_segmentation, output_prefix, padding=10, labels_to_register=[2, 3, 4, 5], total_sigma=0.5, is_test=False) 2875 2876 if not is_test: 2877 tbinseg = ants.mask_image( cit168labT, cit168labT, [7,9,23,25,33,34], binarize=True) 2878 tbinseg = ants.morphology( tbinseg, "dilate", 14 ) 2879 ibinseg = ants.apply_transforms( img, tbinseg, cit168reg['invtransforms'], 2880 interpolator='nearestNeighbor') 2881 snreg = region_reg( img, myparc['tissue_segmentation'], ibinseg, 2882 cit168adni, tbinseg, output_prefix=output_prefix + "_SNREG", 2883 padding = 4, is_test=False )['synL'] 2884 tbinseg = ants.mask_image( cit168labT, cit168labT, [7,9,23,25,33,34], binarize=False) 2885 snseg = ants.apply_transforms( img, tbinseg, 2886 snreg['invtransforms'], interpolator = 'nearestNeighbor' ) 2887 snseg = snseg * ants.threshold_image( myparc['tissue_segmentation'], 2, 6 ) 2888 snseg_desc = map_segmentation_to_dataframe( 'CIT168_Reinf_Learn_v1_label_descriptions_pad', snseg ).dropna(axis=0) 2889 else : 2890 snseg = None 2891 snseg_desc = None 2892 2893 if verbose: 2894 print( "brainstem and cerebellar segmentation" ) 2895 # midbrain/brainstem 2896 brainstemseg = ants.apply_transforms( img, cit168labStem, 2897 cit168reg['invtransforms'], interpolator = 'nearestNeighbor' ) 2898 brainstemseg = brainstemseg * braintissuemask 2899 brainstem_desc = map_segmentation_to_dataframe( 'CIT168_T1w_700um_pad_adni_brainstem', brainstemseg ) 2900 brainstem_desc = brainstem_desc.loc[:, ~brainstem_desc.columns.str.contains('^Side')] 2901 2902 2903 # cerebellum 2904 cereb = antspynet.cerebellum_morphology( ants.iMath( x, "Normalize" ), compute_thickness_image=False, verbose=False, do_preprocessing=True ) 2905 # refinement with cerebellum estimate (comment out since it's not needed for this image). 2906 maskc = ants.threshold_image(cereb['cerebellum_probability_image'], 0.5, 1, 1, 0) 2907 cereb = antspynet.cerebellum_morphology( ants.iMath( x, "Normalize" ), cerebellum_mask=maskc, compute_thickness_image=False, verbose=False, do_preprocessing=True ) 2908 cereb_desc = map_segmentation_to_dataframe( 'cerebellum', cereb['parcellation_segmentation_image'] ).dropna(axis=0) 2909 2910 if verbose: 2911 print( "antspyt1w.hierarchical complete" ) 2912 mydataframes = { 2913 "rbp": myqc['brain'], 2914 "icv": myicv, 2915 "hemispheres":hemi, 2916 "tissues":tissue, 2917 "dktlobes":dktl, 2918 "dktregions":dktp, 2919 "dktcortex":dktc, 2920 "wmtracts_left":wmtdfL, 2921 "wmtracts_right":wmtdfR, 2922# "wmh":myhypo['wmh_summary'], 2923 "mtl":deep_flash['mtl_description'], 2924 "bf":deep_bf['description'], 2925 "cit168":cit168lab_desc, 2926 "deep_cit168":deep_cit['description'], 2927 "snseg":snseg_desc, 2928# "hippLR":hippLR['description'], 2929 "brainstem": brainstem_desc, 2930 "cerebellum": cereb_desc 2931 } 2932 2933 outputs = { 2934 "brain_n4_dnz": img, 2935 "brain_n4_dnz_png": myqc['brain_image'], 2936 "brain_extraction": imgbxt, 2937 "tissue_seg_png": tissue_seg_png, 2938 "left_right": mylr, 2939 "dkt_parc": myparc, 2940 "registration":reg, 2941# "hippLR":hippLR['segmentation'], 2942# "white_matter_hypointensity":myhypo['wmh_probability_image'], 2943 "wm_tractsL":wm_tractsL, 2944 "wm_tractsR":wm_tractsR, 2945 "mtl":deep_flash['mtl_segmentation'], 2946 "bf":deep_bf['segmentation'], 2947 "deep_cit168lab": deep_cit['segmentation'], 2948 "cit168lab": cit168lab, 2949 "cit168reg": cit168reg, 2950 "snseg":snseg, 2951 "snreg":snreg, 2952 "brainstem": brainstemseg, 2953 "cerebellum":cereb['parcellation_segmentation_image'], 2954 "dataframes": mydataframes 2955 } 2956 2957 return outputs 2958 2959 2960def trim_segmentation_by_distance( segmentation, which_label, distance ): 2961 """ 2962 trim a segmentation by the distance provided by the user. computes a distance 2963 transform from the segmentation - treated as binary - and trims the target 2964 label by that distance. 2965 2966 Arguments 2967 --------- 2968 segmentation : ants image segmentation 2969 2970 which_label : the label to trim 2971 2972 distance : float distance value 2973 2974 Returns 2975 ------- 2976 trimmed_segmentation 2977 2978 Example 2979 ------- 2980 >>> import ants 2981 >>> img = ants.image_read( ants.get_data( 'r16' ) ) 2982 >>> seg = ants.threshold_image( img, "Otsu", 3 ) 2983 >>> tseg = antspyt1w.trim_segmentation_by_distance( seg, 1, 10 ) 2984 """ 2985 bseg = ants.threshold_image( segmentation, 1, segmentation.max() ) 2986 dist = ants.iMath( bseg, "MaurerDistance" ) * (-1.0) 2987 disttrim = ants.threshold_image( dist, distance, dist.max() ) 2988 tarseg = ants.threshold_image( segmentation, which_label, which_label ) * disttrim 2989 segmentationtrim = segmentation.clone() 2990 segmentationtrim[ segmentation == which_label ] = 0 2991 return segmentationtrim + tarseg * which_label 2992 2993 2994 2995def zoom_syn( target_image, template, template_segmentations, 2996 initial_registration, 2997 dilation = 4, 2998 regIterations = [25] ): 2999 """ 3000 zoomed in syn - a hierarchical registration applied to a hierarchical segmentation 3001 3002 Initial registration is followed up by a refined and focused high-resolution registration. 3003 This is performed on the cropped image where the cropping region is determined 3004 by the first segmentation in the template_segmentations list. Segmentations 3005 after the first one are assumed to exist as sub-regions of the first. All 3006 segmentations are assumed to be binary. 3007 3008 Arguments 3009 --------- 3010 target_image : ants image at original resolution 3011 3012 template : ants image template to be mapped to the target image 3013 3014 template_segmentations : list of binary segmentation images 3015 3016 dilation : morphological dilation amount applied to the first segmentation and used for cropping 3017 3018 regIterations : parameter passed to ants.registration 3019 3020 Returns 3021 ------- 3022 dictionary 3023 containing segmentation and registration results in addition to cropping results 3024 3025 Example 3026 ------- 3027 >>> import ants 3028 >>> ireg = ants.registration( target_image, template, "antsRegistrationSyNQuickRepro[s]" ) 3029 >>> xxx = antspyt1w.zoom_syn( orb, template, level2segs, ireg ) 3030 """ 3031 croppertem = ants.iMath( template_segmentations[0], "MD", dilation ) 3032 templatecrop = ants.crop_image( template, croppertem ) 3033 cropper = ants.apply_transforms( target_image, 3034 croppertem, initial_registration['fwdtransforms'], 3035 interpolator='linear' ).threshold_image(0.5,1.e9) 3036 croplow = ants.crop_image( target_image, cropper ) 3037 synnerlow = ants.registration( croplow, templatecrop, 3038 'SyNOnly', gradStep = 0.20, regIterations = regIterations, randomSeed=1, 3039 syn_metric='cc', syn_sampling=2, total_sigma=0.5, 3040 initialTransform = initial_registration['fwdtransforms'] ) 3041 orlist = [] 3042 for jj in range(len(template_segmentations)): 3043 target_imageg = ants.apply_transforms( target_image, template_segmentations[jj], 3044 synnerlow['fwdtransforms'], 3045 interpolator='linear' ).threshold_image(0.5,1e9) 3046 orlist.append( target_imageg ) 3047 return{ 3048 'segmentations': orlist, 3049 'registration': synnerlow, 3050 'croppedimage': croplow, 3051 'croppingmask': cropper 3052 } 3053 3054 3055def read_hierarchical( output_prefix ): 3056 """ 3057 standardized reading of output for hierarchical function 3058 3059 Arguments 3060 --------- 3061 3062 output_prefix : string path including directory and file prefix that was 3063 be applied to all output, both csv and images. 3064 3065 Returns 3066 ------- 3067 hierarchical_object : output of antspyt1w.hierarchical 3068 objects will be None if they cannot be found or if reading/writing 3069 is not yet implemented for them 3070 3071 """ 3072 3073 mydataframes = { 3074 "rbp": None, 3075 "hemispheres":None, 3076 "tissues":None, 3077 "dktlobes":None, 3078 "dktregions":None, 3079 "dktcortex":None, 3080 "wmtracts_left":None, 3081 "wmtracts_right":None, 3082 # "wmh":myhypo['wmh_summary'], 3083 "mtl":None, 3084 "bf":None, 3085 "cit168":None, 3086 "deep_cit168":None, 3087 "snseg":None, 3088# "hippLR":None, 3089 "cerebellum":None, 3090 "brainstem":None 3091 } 3092 3093 dkt_parc = { 3094 "tissue_segmentation":None, 3095 "tissue_probabilities":None, 3096 "dkt_parcellation":None, 3097 "dkt_lobes":None, 3098 "dkt_cortex": None, 3099 "hemisphere_labels": None, 3100 "wmSNR": None, 3101 "wmcsfSNR": None, } 3102 3103 hierarchical_object = { 3104 "brain_n4_dnz": None, 3105 "brain_n4_dnz_png": None, 3106 "brain_extraction": None, 3107 "tissue_seg_png": None, 3108 "left_right": None, 3109 "dkt_parc": dkt_parc, 3110 "registration":None, 3111# "hippLR":None, 3112 # "white_matter_hypointensity":None, 3113 "wm_tractsL":None, 3114 "wm_tractsR":None, 3115 "mtl":None, 3116 "bf":None, 3117 "deep_cit168lab": None, 3118 "cit168lab": None, 3119 "cit168reg": None, 3120 "snseg":None, 3121 "snreg":None, 3122 "cerebellum":None, 3123 "brainstem":None, 3124 "dataframes": mydataframes 3125 } 3126 3127 for myvar in hierarchical_object['dataframes'].keys(): 3128 if hierarchical_object['dataframes'][myvar] is None and exists( output_prefix + myvar + ".csv"): 3129 hierarchical_object['dataframes'][myvar] = pd.read_csv(output_prefix + myvar + ".csv") 3130 3131 myvarlist = hierarchical_object.keys() 3132 for myvar in myvarlist: 3133 if hierarchical_object[myvar] is None and exists( output_prefix + myvar + '.nii.gz' ): 3134 hierarchical_object[myvar] = ants.image_read( output_prefix + myvar + '.nii.gz' ) 3135 3136 myvarlist = [ 3137 'tissue_segmentation', 3138 'dkt_parcellation', 3139 'dkt_lobes', 3140 'dkt_cortex', 3141 'hemisphere_labels' ] 3142 for myvar in myvarlist: 3143 if hierarchical_object['dkt_parc'][myvar] is None and exists( output_prefix + myvar + '.nii.gz' ): 3144 hierarchical_object['dkt_parc'][myvar] = ants.image_read( output_prefix + myvar + '.nii.gz' ) 3145 3146 3147 return hierarchical_object 3148 3149 3150 3151 3152def write_hierarchical( hierarchical_object, output_prefix, verbose=False ): 3153 """ 3154 standardized writing of output for hierarchical function 3155 3156 Arguments 3157 --------- 3158 hierarchical_object : output of antspyt1w.hierarchical 3159 3160 output_prefix : string path including directory and file prefix that will 3161 be applied to all output, both csv and images. 3162 3163 verbose: boolean 3164 3165 Returns 3166 ------- 3167 None 3168 3169 """ 3170 3171 # write extant dataframes 3172 for myvar in hierarchical_object['dataframes'].keys(): 3173 if hierarchical_object['dataframes'][myvar] is not None: 3174 if verbose: 3175 print( myvar ) 3176 hierarchical_object['dataframes'][myvar].dropna(axis=0).to_csv(output_prefix + myvar + ".csv") 3177 3178 myvarlist = hierarchical_object.keys() 3179 r16img = ants.image_read( ants.get_data( "r16" )) 3180 for myvar in myvarlist: 3181 if hierarchical_object[myvar] is not None and type(hierarchical_object[myvar]) == type( r16img ): 3182 if verbose: 3183 print( output_prefix + myvar ) 3184 ants.image_write( hierarchical_object[myvar], output_prefix + myvar + '.nii.gz' ) 3185 3186 myvarlist = [ 3187 'tissue_segmentation', 3188 'dkt_parcellation', 3189 'dkt_lobes', 3190 'dkt_cortex', 3191 'hemisphere_labels' ] 3192 for myvar in myvarlist: 3193 if hierarchical_object['dkt_parc'][myvar] is not None: 3194 if verbose: 3195 print( output_prefix + myvar ) 3196 ants.image_write( hierarchical_object['dkt_parc'][myvar], output_prefix + myvar + '.nii.gz' ) 3197 3198 return 3199 3200 3201 3202 3203def merge_hierarchical_csvs_to_wide_format( hierarchical_dataframes, col_names = None , identifier=None, identifier_name='u_hier_id', verbose=False ): 3204 """ 3205 standardized merging of output for dataframes produced by hierarchical function. 3206 3207 Arguments 3208 --------- 3209 hierarchical_dataframes : output of antspyt1w.hierarchical 3210 3211 identifier : unique subject identifier e.g. subject_001 3212 3213 identifier_name : string name for the unique identifer column e.g. subject_id 3214 3215 Returns 3216 ------- 3217 data frame in wide format 3218 3219 """ 3220 if identifier is None: 3221 identifier='A' 3222 wide_df = pd.DataFrame( ) 3223 icvkey='icv' 3224 if icvkey in hierarchical_dataframes.keys(): 3225 temp = hierarchical_dataframes[icvkey].copy() 3226 temp = temp.loc[:, ~temp.columns.str.contains('^Unnamed')] 3227 temp.insert(loc=0, column=identifier_name, value=identifier) 3228 temp = temp.set_index(identifier_name) 3229 wide_df = wide_df.join(temp,how='outer') 3230 for myvar in hierarchical_dataframes.keys(): 3231 if verbose: 3232 print( myvar ) 3233 if hierarchical_dataframes[myvar] is not None: 3234 jdf = hierarchical_dataframes[myvar].dropna(axis=0) 3235 jdf = jdf.loc[:, ~jdf.columns.str.contains('^Unnamed')] 3236 if col_names is not None : 3237 for col_name in col_names : 3238 if jdf.shape[0] > 1 and any( jdf.columns.str.contains(col_name)): 3239 varsofinterest = ["Description", col_name] 3240 jdfsub = jdf[varsofinterest] 3241 jdfsub.insert(loc=0, column=identifier_name, value=identifier) 3242 jdfsub = jdfsub.set_index([identifier_name, 'Description'])[col_name].unstack().add_prefix(col_name + '_') 3243 jdfsub.columns=jdfsub.columns 3244 jdfsub = jdfsub.rename(mapper=lambda x: x.strip().replace(' ', '_').lower(), axis=1) 3245 wide_df = wide_df.join(jdfsub,how='outer') 3246 if jdf.shape[0] > 1 and any( jdf.columns.str.contains('VolumeInMillimeters')): 3247 varsofinterest = ["Description", "VolumeInMillimeters"] 3248 jdfsub = jdf[varsofinterest] 3249 jdfsub.insert(loc=0, column=identifier_name, value=identifier) 3250 jdfsub=jdfsub.set_index([identifier_name, 'Description']).VolumeInMillimeters.unstack().add_prefix('Vol_') 3251 jdfsub.columns=jdfsub.columns+myvar 3252 jdfsub = jdfsub.rename(mapper=lambda x: x.strip().replace(' ', '_').lower(), axis=1) 3253 wide_df = wide_df.join(jdfsub,how='outer') 3254 if jdf.shape[0] > 1 and any( jdf.columns.str.contains('SurfaceAreaInMillimetersSquared')): 3255 varsofinterest = ["Description", "SurfaceAreaInMillimetersSquared"] 3256 jdfsub = jdf[varsofinterest] 3257 jdfsub.insert(loc=0, column=identifier_name, value=identifier) 3258 jdfsub=jdfsub.set_index([identifier_name, 'Description']).SurfaceAreaInMillimetersSquared.unstack().add_prefix('Area_') 3259 jdfsub.columns=jdfsub.columns+myvar 3260 jdfsub = jdfsub.rename(mapper=lambda x: x.strip().replace(' ', '_').lower(), axis=1) 3261 wide_df = wide_df.join(jdfsub,how='outer') 3262 if jdf.shape[0] > 1 and any( jdf.columns.str.contains('SurfaceAreaInMillimetersSquared')) and any( jdf.columns.str.contains('VolumeInMillimeters')): 3263 varsofinterest = ["Description", "VolumeInMillimeters", "SurfaceAreaInMillimetersSquared"] 3264 jdfsub = jdf[varsofinterest] 3265 jdfsub.insert(loc=0, column=identifier_name, value=identifier) 3266 jdfsub.insert(loc=1, column='thickness',value=jdfsub['VolumeInMillimeters']/jdfsub['SurfaceAreaInMillimetersSquared']) 3267 jdfsub=jdfsub.set_index([identifier_name, 'Description']).thickness.unstack().add_prefix('Thk_') 3268 jdfsub.columns=jdfsub.columns+myvar 3269 jdfsub = jdfsub.rename(mapper=lambda x: x.strip().replace(' ', '_').lower(), axis=1) 3270 wide_df = wide_df.join(jdfsub,how='outer') 3271 rbpkey='rbp' 3272 if rbpkey in hierarchical_dataframes.keys(): 3273 temp = hierarchical_dataframes[rbpkey].copy() 3274 temp = temp.loc[:, ~temp.columns.str.contains('^Unnamed')] 3275 temp.insert(loc=0, column=identifier_name, value=identifier) 3276 temp = temp.set_index(identifier_name) 3277 wide_df = wide_df.join(temp,how='outer') 3278 # handle wmh 3279 wmhkey='wmh' 3280 if wmhkey in hierarchical_dataframes.keys(): 3281 df=hierarchical_dataframes[wmhkey].copy() 3282 df.insert(loc=0, column=identifier_name, value=identifier) 3283 df = df.set_index(identifier_name) 3284 wmhval = df.loc[ df.Description == 'Volume_of_WMH','Value'] 3285 wide_df.insert(loc = 0, column = 'wmh_vol', value =wmhval ) 3286 wmhval = df.loc[ df.Description == 'Integral_WMH_probability','Value'] 3287 wide_df.insert(loc = 0, column = 'wmh_integral_prob', value =wmhval ) 3288 wmhval = df.loc[ df.Description == 'Log_Evidence','Value'] 3289 wide_df.insert(loc = 0, column = 'wmh_log_evidence', value =wmhval ) 3290 wide_df['wmh_log_evidence']=wmhval 3291 wide_df.insert(loc = 0, column = identifier_name, value = identifier) 3292 return wide_df 3293 3294 3295 3296 3297def super_resolution_segmentation_per_label( 3298 imgIn, 3299 segmentation, # segmentation label image 3300 upFactor, 3301 sr_model, 3302 segmentation_numbers, 3303 dilation_amount = 0, 3304 probability_images=None, # probability list 3305 probability_labels=None, # the segmentation ids for the probability image, 3306 max_lab_plus_one=True, 3307 target_range=[1,0], 3308 poly_order = 'hist', 3309 verbose = False, 3310): 3311 """ 3312 Apply a two-channel super resolution model to an image and probability pair. 3313 3314 Arguments 3315 --------- 3316 imgIn : ANTsImage 3317 image to be upsampled 3318 3319 segmentation : ANTsImage 3320 segmentation probability, n-ary or binary image 3321 3322 upFactor : list 3323 the upsampling factors associated with the super-resolution model 3324 3325 sr_model : tensorflow model or String 3326 for computing super-resolution; String denotes an interpolation type 3327 3328 segmentation_numbers : list of target segmentation labels 3329 list containing integer segmentation labels 3330 3331 dilation_amount : integer 3332 amount to pad around each segmentation defining region over which to 3333 computer the super-resolution in each label 3334 3335 probability_images : list of ANTsImages 3336 list of probability images 3337 3338 probability_labels : integer list 3339 providing the integer values associated with each probability image 3340 3341 max_lab_plus_one : boolean 3342 add background label 3343 3344 target_range : lower and upper limit of intensity expected by network [1,0] if trained by siq 3345 3346 poly_order : integer or 'hist' or None 3347 3348 verbose : boolean 3349 whether to show status updates 3350 3351 Returns 3352 ------- 3353 3354 dictionary w/ following key/value pairs: 3355 `super_resolution` : ANTsImage 3356 super_resolution image 3357 3358 `super_resolution_segmentation` : ANTsImage 3359 super_resolution_segmentation image 3360 3361 `segmentation_geometry` : list of data frame types 3362 segmentation geometry for each label 3363 3364 `probability_images` : list of ANTsImage 3365 segmentation probability maps 3366 3367 3368 Example 3369 ------- 3370 >>> import ants 3371 >>> ref = ants.image_read( ants.get_ants_data('r16')) 3372 """ 3373 import re 3374 if type( sr_model ) == type(""): 3375 if re.search( 'h5', sr_model ) is not None: 3376 if verbose: 3377 print("load model") 3378 sr_model=tf.keras.models.load_model( sr_model, compile=False ) 3379 newspc = ( np.asarray( ants.get_spacing( imgIn ) ) ).tolist() 3380 for k in range(len(newspc)): 3381 newspc[k] = newspc[k]/upFactor[k] 3382 imgup = ants.resample_image( imgIn, newspc, use_voxels=False, interp_type=0 ) 3383 imgsrfull = imgup * 0.0 3384 weightedavg = imgup * 0.0 3385 problist = [] 3386 bkgdilate = 2 3387 segmentationUse = ants.image_clone( segmentation ) 3388 segmentationUse = ants.mask_image( segmentationUse, segmentationUse, segmentation_numbers ) 3389 segmentation_numbers_use = segmentation_numbers.copy() 3390 if max_lab_plus_one: 3391 background = ants.mask_image( segmentationUse, segmentationUse, segmentation_numbers, binarize=True ) 3392 background = ants.iMath(background,"MD",bkgdilate) - background 3393 backgroundup = ants.resample_image_to_target( background, imgup, interp_type='linear' ) 3394 segmentation_numbers_use.append( max(segmentation_numbers) + 1 ) 3395 segmentationUse = segmentationUse + background * max(segmentation_numbers_use) 3396 3397 for locallab in segmentation_numbers: 3398 if verbose: 3399 print( "SR-per-label:" + str( locallab ) ) 3400 binseg = ants.threshold_image( segmentationUse, locallab, locallab ) 3401 sizethresh = 2 3402 if ( binseg == 1 ).sum() < sizethresh : 3403 warnings.warn( "SR-per-label:" + str( locallab ) + ' is small' ) 3404 if ( binseg == 1 ).sum() >= sizethresh : 3405 minprob="NA" 3406 maxprob="NA" 3407 if probability_images is not None: 3408 whichprob = probability_labels.index(locallab) 3409 probimg = probability_images[whichprob].resample_image_to_target( binseg ) 3410 minprob = min( probimg[ binseg >= 0.5 ] ) 3411 maxprob = max( probimg[ binseg >= 0.5 ] ) 3412 if verbose: 3413 print( "SR-per-label:" + str( locallab ) + " min/max-prob: " + str(minprob)+ " / " + str(maxprob) ) 3414 binsegdil = ants.iMath( ants.threshold_image( segmentationUse, locallab, locallab ), "MD", dilation_amount ) 3415 binsegdil2input = ants.resample_image_to_target( binsegdil, imgIn, interp_type='nearestNeighbor' ) 3416 imgc = ants.crop_image( ants.iMath(imgIn,"Normalize"), binsegdil2input ) 3417 imgc = imgc * target_range[0] - target_range[1] # for SR 3418 imgchCore = ants.crop_image( binseg, binsegdil ) 3419 if probability_images is not None: 3420 imgchCore = ants.crop_image( probimg, binsegdil ) 3421# imgch = ants.iMath( imgchCore, "Normalize" ) * target_range[0] - target_range[1] # for SR 3422 imgch = imgchCore * target_range[0] - target_range[1] # for SR 3423 if type( sr_model ) == type(""): # this is just for testing 3424 binsegup = ants.resample_image_to_target( binseg, imgup, interp_type='linear' ) 3425 problist.append( binsegup ) 3426 else: 3427 if verbose: 3428 print(imgc) 3429 print(imgch) 3430 myarr = np.stack( [imgc.numpy(),imgch.numpy()],axis=3 ) 3431 newshape = np.concatenate( [ [1],np.asarray( myarr.shape )] ) 3432 myarr = myarr.reshape( newshape ) 3433 if verbose: 3434 print("calling prediction function") 3435 print("arr range " + str(myarr.max() ) + " " + str(myarr.min() ) ) 3436 print( myarr.shape ) 3437 pred = sr_model.predict( myarr ) 3438 predshape = np.asarray(pred).shape 3439 predshapelen = len( predshape ) 3440 if predshape[ predshapelen - 1] == 2: 3441 pred = tf.split( pred, 2, axis=4 ) 3442 if verbose: 3443 print("predict done") 3444 imgsr = ants.from_numpy( tf.squeeze( pred[0] ).numpy()) 3445 if verbose: 3446 print( imgc ) 3447 print( imgsr ) 3448 imgsr = ants.copy_image_info( imgc, imgsr ) 3449 ants.set_spacing( imgsr, newspc ) 3450 imgsrh = ants.from_numpy( tf.squeeze( pred[1] ).numpy()) 3451 imgsrh = ants.copy_image_info( imgc, imgsrh ) 3452 ants.set_spacing( imgsrh, newspc ) 3453 if poly_order is not None: 3454 if poly_order == 'hist': 3455 if verbose: 3456 print("hist match intensity") 3457 imgsrh = ants.histogram_match_image( imgsrh, imgchCore ) 3458 imgsr = ants.histogram_match_image( imgsr, imgc ) 3459 else: 3460 if verbose: 3461 print("poly match intensity") 3462 imgsr = ants.regression_match_image( imgsr, ants.resample_image_to_target(imgup,imgsr), poly_order=poly_order ) 3463 imgsrh = ants.regression_match_image( imgsrh, ants.resample_image_to_target(imgchCore,imgsr), poly_order=poly_order ) 3464 problist.append( imgsrh ) 3465 contribtoavg = ants.resample_image_to_target( imgsr*0+1, imgup, interp_type='nearestNeighbor' ) 3466 weightedavg = weightedavg + contribtoavg 3467 imgsrfull = imgsrfull + ants.resample_image_to_target( imgsr, imgup, interp_type='nearestNeighbor' ) 3468 3469 if max_lab_plus_one: 3470 problist.append( backgroundup ) 3471 3472 imgsrfull2 = imgsrfull 3473 selector = imgsrfull == 0 3474 imgsrfull2[ selector ] = imgup[ selector ] 3475 weightedavg[ weightedavg == 0.0 ] = 1.0 3476 imgsrfull2=imgsrfull2/weightedavg 3477 imgsrfull2[ imgup == 0 ] = 0 3478 3479 for k in range(len(problist)): 3480 problist[k] = ants.resample_image_to_target(problist[k],imgsrfull2,interp_type='linear') 3481 3482 if max_lab_plus_one: 3483 tarmask = ants.threshold_image( segmentationUse, 1, segmentationUse.max() ) 3484 else: 3485 tarmask = ants.threshold_image( segmentationUse, 1, segmentationUse.max() ).iMath("MD",1) 3486 tarmask = ants.resample_image_to_target( tarmask, imgsrfull2, interp_type='nearestNeighbor' ) 3487 3488 3489 if True: # we dont generally know which of these is background/foreground so ..... 3490 segmat = ants.images_to_matrix(problist, tarmask) 3491 finalsegvec = segmat.argmax(axis=0) 3492 finalsegvec2 = finalsegvec.copy() 3493 for i in range(len(problist)): 3494 segnum = segmentation_numbers_use[i] 3495 finalsegvec2[finalsegvec == i] = segnum 3496 outimg = ants.make_image(tarmask, finalsegvec2) 3497 outimg = ants.mask_image( outimg, outimg, segmentation_numbers ) 3498 seggeom = ants.label_geometry_measures( outimg ) 3499 else: 3500 image_matrix = ants.image_list_to_matrix(problist[0:(len(problist)-1)], tarmask ) 3501 background_foreground_matrix = np.stack([ants.image_list_to_matrix([problist[len(problist)-1]], tarmask), 3502 np.expand_dims(np.sum(image_matrix, axis=0), axis=0)]) 3503 foreground_matrix = np.argmax(background_foreground_matrix, axis=0) 3504 segmentation_matrix = (np.argmax(image_matrix, axis=0) + 1) * foreground_matrix 3505 segmentation_image = ants.matrix_to_images(np.expand_dims(segmentation_matrix, axis=0), tarmask)[0] 3506 outimg = ants.image_clone(segmentation_image) 3507 for i in range(len(problist)): 3508 outimg[segmentation_image==i] = segmentation_numbers_use[i] 3509 outimg = ants.mask_image( outimg, outimg, segmentation_numbers ) 3510 seggeom = ants.label_geometry_measures( outimg ) 3511 3512 return { 3513 "super_resolution": imgsrfull2, 3514 "super_resolution_segmentation": outimg, 3515 "segmentation_geometry": seggeom, 3516 "probability_images": problist 3517 } 3518 3519 3520 3521 3522def super_resolution_segmentation_with_probabilities( 3523 img, 3524 initial_probabilities, 3525 sr_model, 3526 target_range=[1,0], 3527 verbose = False 3528): 3529 """ 3530 Simultaneous super-resolution and probabilistic segmentation. 3531 3532 this function will perform localized super-resolution analysis on the 3533 underlying image content to produce a SR version of both the local region 3534 intentsity and segmentation probability. the number of probability images 3535 determines the number of outputs for the intensity and probabilitiy list 3536 that are returned. NOTE: should consider inverse logit transform option. 3537 3538 Arguments 3539 --------- 3540 img : input intensity image 3541 3542 initial_probabilities: original resolution probability images 3543 3544 sr_model : the super resolution model - should have 2 channel input 3545 3546 target_range : lower and upper limit of intensity expected by network [1,0] if trained by siq 3547 3548 verbose : boolean 3549 3550 Returns 3551 ------- 3552 list of ANTsImages 3553 3554 """ 3555 3556 # SR Part 3557 srimglist = [] 3558 srproblist = [] 3559 mypt = 1.0 / len(initial_probabilities) 3560 3561 for k in range(len(initial_probabilities)): 3562 if k == 0: 3563 srimglist.append( img ) 3564 srproblist.append( initial_probabilities[k] ) 3565 else: 3566 tempm = ants.threshold_image( initial_probabilities[k], mypt, 2 ).iMath("MD",1) 3567 imgc = ants.crop_image(img,tempm) 3568 imgch = ants.crop_image(initial_probabilities[k],tempm) 3569 if verbose: 3570 print(k) 3571 print(imgc) 3572 imgcrescale = ants.iMath( imgc, "Normalize" ) * target_range[0] - target_range[1] # for SR 3573 imgchrescale = imgch * target_range[0] - target_range[1] 3574 myarr = np.stack( [ imgcrescale.numpy(), imgchrescale.numpy() ],axis=3 ) 3575 newshape = np.concatenate( [ [1],np.asarray( myarr.shape )] ) 3576 myarr = myarr.reshape( newshape ) 3577 pred = sr_model.predict( myarr ) 3578 imgsr = ants.from_numpy( tf.squeeze( pred[0] ).numpy()) 3579 imgsr = ants.copy_image_info( imgc, imgsr ) 3580 newspc = ( np.asarray( ants.get_spacing( imgsr ) ) * 0.5 ).tolist() 3581 ants.set_spacing( imgsr, newspc ) 3582 if verbose: 3583 print(imgsr) 3584 imgsr = ants.regression_match_image( imgsr, ants.resample_image_to_target(imgc,imgsr) ) 3585 imgsrh = ants.from_numpy( tf.squeeze( pred[1] ).numpy()) 3586 imgsrh = ants.copy_image_info( imgsr, imgsrh ) 3587 tempup = ants.resample_image_to_target( tempm, imgsr ) 3588 srimglist.append( imgsr ) 3589 # NOTE: get rid of pixellated junk/artifacts - acts like a prior 3590 srproblist.append( imgsrh * tempup ) 3591 3592 if verbose: 3593 print("done srwithprob") 3594 3595 labels = { 3596 'sr_intensities':srimglist, 3597 'sr_probabilities':srproblist, 3598 } 3599 return labels 3600 3601def kelly_kapowski_thickness( x, labels, 3602 label_description='dkt', iterations=45, max_thickness=6.0, mydap=None, verbose=False ): 3603 """ 3604 Apply a two-channel super resolution model to an image and probability pair. 3605 3606 Arguments 3607 --------- 3608 x : ANTsImage 3609 raw t1 image 3610 3611 labels : ANTsImage 3612 cortical parcellation 3613 3614 label_description : pandas data frame 3615 3616 iterations : integer 3617 number of iterations ( probably not to be changed except for testing ) 3618 3619 max_thickness : in mm to consider default is 6.0 3620 3621 verbose : boolean 3622 3623 Returns 3624 ------- 3625 3626 dictionary w/ following key/value pairs: 3627 `thickness_image` : ANTsImage 3628 3629 `thickness_dataframe` : dataframe with mean thickness values 3630 3631 Example 3632 ------- 3633 >>> import ants 3634 >>> thk = antspyt1w.kelly_kapowski_thickness( x ) 3635 """ 3636 if verbose: 3637 myverb=1 3638 else: 3639 myverb=0 3640 seg = deep_tissue_segmentation( x ) 3641 kkthk = ants.kelly_kapowski( s=seg['segmentation_image'], 3642 g=seg['probability_images'][2], w=seg['probability_images'][3], 3643 its=iterations, r=0.025, m=1.5, verbose=myverb ) 3644 kkthkmask = ants.threshold_image( kkthk, 0.25, max_thickness ) 3645 kkdf = map_intensity_to_dataframe( 3646 label_description, 3647 kkthk, 3648 labels * kkthkmask ) 3649 kkdf_wide = merge_hierarchical_csvs_to_wide_format( {'KK' : kkdf}, col_names = ['Mean'] ) 3650 return { 3651 'thickness_image' : kkthk, 3652 'thickness_dataframe' : kkdf_wide 3653 } 3654 3655 3656def minimal_sr_preprocessing( x, imgbxt=None ): 3657 """ 3658 x : input t1 image (raw) 3659 3660 imgbxt : optional existing brain extraction 3661 3662 outputs: preprocessedT1, hemisphereLabels 3663 """ 3664 if x.dimension != 3: 3665 raise ValueError('hierarchical: input image should be 3-dimensional') 3666 3667 tfn = get_data('T_template0', target_extension='.nii.gz' ) 3668 tlrfn = get_data('T_template0_LR', target_extension='.nii.gz' ) 3669 3670 ##### read images and do simple bxt ops 3671 templatea = ants.image_read( tfn ) 3672 templatea = ( templatea * antspynet.brain_extraction( templatea, 't1' ) ).iMath( "Normalize" ) 3673 templatealr = ants.image_read( tlrfn ) 3674 if imgbxt is None: 3675 imgbxt = brain_extraction( ants.iMath( x, "Normalize" ) ) 3676 img = preprocess_intensity( ants.iMath( x, "Normalize" ), imgbxt ) 3677 img = ants.iMath( img, "Normalize" ) 3678 mylr = label_hemispheres( img, templatea, templatealr ) 3679 return img, mylr * imgbxt
def
get_data(name=None, force_download=False, version=46, target_extension='.csv'):
112def get_data(name=None, force_download=False, version=46, target_extension='.csv'): 113 """ 114 Get ANTsPyT1w data filename 115 116 The first time this is called, it will download data to ~/.antspyt1w. 117 After, it will just read data from disk. The ~/.antspyt1w may need to 118 be periodically deleted in order to ensure data is current. 119 120 Arguments 121 --------- 122 name : string 123 name of data tag to retrieve 124 Options: 125 - 'all' 126 - 'dkt' 127 - 'hemisphere' 128 - 'lobes' 129 - 'tissues' 130 - 'T_template0' 131 - 'T_template0_LR' 132 - 'T_template0_LobesBstem' 133 - 'T_template0_WMP' 134 - 'T_template0_Symmetric' 135 - 'T_template0_SymmetricLR' 136 - 'PPMI-3803-20120814-MRI_T1-I340756' 137 - 'simwmseg' 138 - 'simwmdisc' 139 - 'wmh_evidence' 140 - 'wm_major_tracts' 141 142 force_download: boolean 143 144 version: version of data to download (integer) 145 146 Returns 147 ------- 148 string 149 filepath of selected data 150 151 Example 152 ------- 153 >>> import ants 154 >>> ppmi = ants.get_ants_data('ppmi') 155 """ 156 import os 157 import shutil 158 from pathlib import Path 159 import tensorflow as tf 160 DATA_PATH = os.path.join(os.path.expanduser('~'), '.antspyt1w') 161 os.makedirs(DATA_PATH, exist_ok=True) 162 163 def mv_subfolder_files(folder, verbose=False): 164 """ 165 Move files from subfolders to the parent folder. 166 167 Parameters 168 ---------- 169 folder : str 170 Path to the folder. 171 verbose : bool, optional 172 Print information about the files and folders being processed (default is False). 173 174 Returns 175 ------- 176 None 177 """ 178 import os 179 import shutil 180 for root, dirs, files in os.walk(folder): 181 if verbose: 182 print(f"Processing directory: {root}") 183 print(f"Subdirectories: {dirs}") 184 print(f"Files: {files}") 185 186 for file in files: 187 if root != folder: 188 if verbose: 189 print(f"Moving file: {file} from {root} to {folder}") 190 shutil.move(os.path.join(root, file), folder) 191 192 for dir in dirs: 193 if root != folder: 194 if verbose: 195 print(f"Removing directory: {dir} from {root}") 196 shutil.rmtree(os.path.join(root, dir)) 197 198 def download_data(version): 199 url = "https://ndownloader.figshare.com/articles/14766102/versions/" + str(version) 200 target_file_name = "14766102.zip" 201 target_file_name_path = tf.keras.utils.get_file(target_file_name, url, 202 cache_subdir=DATA_PATH, extract=True) 203 os.remove(os.path.join(DATA_PATH, target_file_name)) 204 205 if force_download: 206 download_data(version=version) 207 208 mv_subfolder_files( os.path.expanduser("~/.antspyt1w"), False ) 209 210 # Move files from subdirectories to the main directory 211 for root, dirs, files in os.walk(DATA_PATH): 212 for file in files: 213 if root != DATA_PATH: 214 shutil.move(os.path.join(root, file), DATA_PATH) 215 for dir in dirs: 216 if root != DATA_PATH: 217 shutil.rmtree(os.path.join(root, dir)) 218 219 files = [] 220 for fname in os.listdir(DATA_PATH): 221 if fname.endswith(target_extension): 222 fname = os.path.join(DATA_PATH, fname) 223 files.append(fname) 224 225 if len(files) == 0: 226 download_data(version=version) 227 for fname in os.listdir(DATA_PATH): 228 if fname.endswith(target_extension): 229 fname = os.path.join(DATA_PATH, fname) 230 files.append(fname) 231 232 mv_subfolder_files( os.path.expanduser("~/.antspyt1w"), False ) 233 234 if name == 'all': 235 return files 236 237 datapath = None 238 239 for fname in os.listdir(DATA_PATH): 240 mystem = Path(fname).resolve().stem 241 mystem = Path(mystem).resolve().stem 242 mystem = Path(mystem).resolve().stem 243 if name == mystem and fname.endswith(target_extension): 244 datapath = os.path.join(DATA_PATH, fname) 245 246 if datapath is None: 247 os.listdir(DATA_PATH) 248 warnings.warn("datapath in get_data is None - must be some issue in downloading the data from figshare. try calling get_data in the setup of your system.") 249 return datapath
Get ANTsPyT1w data filename
The first time this is called, it will download data to ~/.antspyt1w. After, it will just read data from disk. The ~/.antspyt1w may need to be periodically deleted in order to ensure data is current.
Arguments
name : string name of data tag to retrieve Options: - 'all' - 'dkt' - 'hemisphere' - 'lobes' - 'tissues' - 'T_template0' - 'T_template0_LR' - 'T_template0_LobesBstem' - 'T_template0_WMP' - 'T_template0_Symmetric' - 'T_template0_SymmetricLR' - 'PPMI-3803-20120814-MRI_T1-I340756' - 'simwmseg' - 'simwmdisc' - 'wmh_evidence' - 'wm_major_tracts'
force_download: boolean
version: version of data to download (integer)
Returns
string filepath of selected data
Example
>>> import ants
>>> ppmi = ants.get_ants_data('ppmi')
def
map_segmentation_to_dataframe(segmentation_type, segmentation_image):
252def map_segmentation_to_dataframe( segmentation_type, segmentation_image ): 253 """ 254 Match the segmentation to its appropriate data frame. We do not check 255 if the segmentation_type and segmentation_image match; this may be indicated 256 by the number of missing values on output eg in column VolumeInMillimeters. 257 258 Arguments 259 --------- 260 segmentation_type : string 261 name of segmentation_type data frame to retrieve 262 Options: 263 - 'dkt' 264 - 'lobes' 265 - 'tissues' 266 - 'hemisphere' 267 - 'wm_major_tracts' 268 269 segmentation_image : antsImage with same values (or mostly the same) as are 270 expected by segmentation_type 271 272 Returns 273 ------- 274 dataframe 275 276 """ 277 mydf_fn = get_data( segmentation_type ) 278 mydf = pd.read_csv( mydf_fn ) 279 mylgo = ants.label_geometry_measures( segmentation_image ) 280 return pd.merge( mydf, mylgo, how='left', on=["Label"] )
Match the segmentation to its appropriate data frame. We do not check if the segmentation_type and segmentation_image match; this may be indicated by the number of missing values on output eg in column VolumeInMillimeters.
Arguments
segmentation_type : string name of segmentation_type data frame to retrieve Options: - 'dkt' - 'lobes' - 'tissues' - 'hemisphere' - 'wm_major_tracts'
segmentation_image : antsImage with same values (or mostly the same) as are expected by segmentation_type
Returns
dataframe
def
hierarchical( x, output_prefix, labels_to_register=[2, 3, 4, 5], imgbxt=None, img6seg=None, cit168=False, is_test=False, atropos_prior=None, sr_model=None, verbose=True):
2611def hierarchical( x, output_prefix, labels_to_register=[2,3,4,5], 2612 imgbxt=None, img6seg=None, cit168 = False, is_test=False, 2613 atropos_prior=None, 2614 sr_model = None, 2615 verbose=True ): 2616 """ 2617 Default processing for a T1-weighted image. See README. 2618 2619 Arguments 2620 --------- 2621 x : T1-weighted neuroimage antsImage 2622 2623 output_prefix: string directory and prefix 2624 2625 labels_to_register: list of integer segmentation labels (of 1 to 6 as defined 2626 by atropos: csf, gm, wm, dgm, brainstem, cerebellum) to define 2627 the tissue types / regions of the brain to register. set to None to 2628 skip registration which will be faster but omit some results. 2629 2630 imgbxt : pre-existing brain extraction - a binary image - will disable some processing 2631 2632 img6seg: pre-existing brain segmentation - a 6-class image ( ANTs standard ) 2633 2634 cit168 : boolean returns labels from CIT168 atlas with high-resolution registration 2635 otherwise, low-resolution regitration is used. 2636 2637 is_test: boolean ( parameters to run more quickly but with low quality ) 2638 2639 atropos_prior: prior weight for atropos post-processing; set to None if you 2640 do not want to use this. will modify the CSF, GM and WM segmentation to 2641 better fit the image intensity at the resolution of the input image. 2642 2643 sr_model: optional will trigger sr-based upsampling in some functions. 2644 2645 verbose: boolean 2646 2647 Returns 2648 ------- 2649 dataframes and associated derived data 2650 2651 - brain_n4_dnz : extracted brain denoised and bias corrected 2652 - brain_extraction : brain mask 2653 - rbp: random basis projection results 2654 - left_right : left righ hemisphere segmentation 2655 - dkt_parc : dictionary object containing segmentation labels 2656 - registration : dictionary object containing registration results 2657 - hippLR : dictionary object containing hippocampus results 2658 - medial_temporal_lobe : dictionary object containing deep_flash (medial temporal lobe parcellation) results 2659 - white_matter_hypointensity : dictionary object containing WMH results 2660 - wm_tractsL : white matter tracts, left 2661 - wm_tractsR : white matter tracts, right 2662 - dataframes : summary data frames 2663 2664 """ 2665 if x.dimension != 3: 2666 raise ValueError('hierarchical: input image should be 3-dimensional') 2667 2668 if verbose: 2669 print("Read") 2670 tfn = get_data('T_template0', target_extension='.nii.gz' ) 2671 tfnw = get_data('T_template0_WMP', target_extension='.nii.gz' ) 2672 tlrfn = get_data('T_template0_LR', target_extension='.nii.gz' ) 2673 bfn = antspynet.get_antsxnet_data( "croppedMni152" ) 2674 2675 ##### read images and do simple bxt ops 2676 templatea = ants.image_read( tfn ) 2677 if verbose: 2678 print("bxt") 2679 templatea = ( templatea * antspynet.brain_extraction( templatea, 't1' ) ).iMath( "Normalize" ) 2680 templateawmprior = ants.image_read( tfnw ) 2681 templatealr = ants.image_read( tlrfn ) 2682 templateb = ants.image_read( bfn ) 2683 templateb = ( templateb * antspynet.brain_extraction( templateb, 't1' ) ).iMath( "Normalize" ) 2684 if imgbxt is None: 2685 probablySR = False 2686 # brain_extraction( ants.iMath( x, "Normalize" ) ) 2687 imgbxt = antspynet.brain_extraction( ants.iMath( x, "Normalize" ), modality="t1threetissue")['segmentation_image'].threshold_image(1,1) 2688 img = preprocess_intensity( ants.iMath( x, "Normalize" ), imgbxt ) 2689 else: 2690 probablySR = True 2691 img = ants.iMath( x, "Normalize" ) 2692 2693 if verbose: 2694 print("QC -- mean of img : "+str(img.mean())) 2695 print(". -- mean of bxt : "+str(imgbxt.mean())) 2696 2697 # this is an unbiased method for identifying predictors that can be used to 2698 # rank / sort data into clusters, some of which may be associated 2699 # with outlierness or low-quality data 2700 templatesmall = ants.resample_image( templateb, (91,109,91), use_voxels=True ) 2701 myqc = inspect_raw_t1( img, output_prefix=output_prefix, option='brain' ) 2702 2703 if verbose: 2704 print("qc results") 2705 print( myqc ) 2706 2707 ##### intensity modifications 2708 img = ants.iMath( img, "Normalize" ) 2709 2710 if verbose: 2711 print("parcellation -- mean of img " + str(img.mean())) 2712 2713 ##### hierarchical labeling 2714 myparc = deep_brain_parcellation( x, templateb, 2715 img6seg = img6seg, 2716 atropos_prior = atropos_prior, 2717 do_cortical_propagation = not is_test, verbose=verbose ) 2718 2719 2720 cit168lab = None 2721 cit168reg = None 2722 cit168lab_desc = None 2723 cit168adni = get_data( "CIT168_T1w_700um_pad_adni",target_extension='.nii.gz') 2724 cit168adni = ants.image_read( cit168adni ).iMath("Normalize") 2725 cit168labT = get_data( "det_atlas_25_pad_LR_adni", target_extension='.nii.gz' ) 2726 cit168labT = ants.image_read( cit168labT ) 2727 cit168labStem = get_data( "CIT168_T1w_700um_pad_adni_brainstem", target_extension='.nii.gz' ) 2728 cit168labStem = ants.image_read( cit168labStem ) 2729 2730 if verbose: 2731 print("cit168-- mean of input image : " + str(img.mean())) 2732 print("cit168-- mean of input imgbxt : " + str(imgbxt.mean())) 2733 print(". -- mean of input segmentation : " + str(myparc['tissue_segmentation'].mean())) 2734 2735 cit168reg = region_reg( 2736 input_image = img, 2737 input_image_tissue_segmentation=myparc['tissue_segmentation'], 2738 input_image_region_segmentation=imgbxt, 2739 input_template=cit168adni, 2740 input_template_region_segmentation=ants.threshold_image( cit168adni, 0.15, 1 ), 2741 output_prefix=output_prefix + "_CIT168RRSYN", 2742 padding=10, 2743 total_sigma=0.5, 2744 is_test=not cit168 )['synL'] 2745 cit168lab = ants.apply_transforms( img, cit168labT, 2746 cit168reg['invtransforms'], interpolator = 'nearestNeighbor' ) 2747 cit168lab_desc = map_segmentation_to_dataframe( 'CIT168_Reinf_Learn_v1_label_descriptions_pad', cit168lab ).dropna(axis=0) 2748 2749 # optional - quick look at result 2750 # bxt_png = output_prefix + "_brain_extraction_dnz_n4_view.png" 2751 # ants.plot(img * 255.0,axis=2,ncol=8,nslices=24, crop=True, black_bg=False, 2752 # filename = bxt_png ) 2753 2754 if verbose: 2755 print("hemi -- mean of input image : " + str(img.mean())) 2756 2757 # assuming data is reasonable quality, we should proceed with the rest ... 2758 mylr = label_hemispheres( img, templatea, templatealr ) 2759 2760 ##### accumulate data into data frames 2761 hemi = map_segmentation_to_dataframe( "hemisphere", myparc['hemisphere_labels'] ) 2762 tissue = map_segmentation_to_dataframe( "tissues", myparc['tissue_segmentation'] ) 2763 dktl = map_segmentation_to_dataframe( "lobes", myparc['dkt_lobes'] ) 2764 dktp = map_segmentation_to_dataframe( "dkt", myparc['dkt_parcellation'] ) 2765 dktc = None 2766 if not is_test: 2767 dktc = map_segmentation_to_dataframe( "dkt", myparc['dkt_cortex'] ) 2768 2769 tissue_seg_png = output_prefix + "_seg.png" 2770 ants.plot( img*255, myparc['tissue_segmentation'], axis=2, nslices=21, ncol=7, 2771 alpha=0.6, filename=tissue_seg_png, 2772 crop=True, black_bg=False ) 2773 2774 myhypo = None 2775# if verbose: 2776# print("WMH") 2777 ##### below here are more exploratory nice to have outputs 2778# myhypo = t1_hypointensity( 2779# img, 2780# myparc['tissue_segmentation'], # segmentation 2781# myparc['tissue_probabilities'][3], # wm posteriors 2782# templatea, 2783# templateawmprior ) 2784 2785 if verbose: 2786 print("icv") 2787 2788 myicv = icv( x ) 2789 2790 braintissuemask = ants.threshold_image( myparc['tissue_segmentation'], 2, 6 ) 2791 if verbose: 2792 print( myicv ) 2793 print("NBM : mean of input image " + str( (img * braintissuemask).mean() ) ) 2794 ##### deep_nbm basal forebrain parcellation 2795 deep_bf = deep_nbm( img * braintissuemask, 2796 get_data("deep_nbm_rank",target_extension='.h5'), 2797 csfquantile=None, aged_template=True, verbose=verbose ) 2798 2799 if is_test: 2800 mydataframes = { 2801 "icv": myicv, 2802 "rbp": myqc['brain'], 2803 "tissues":tissue, 2804 "dktlobes":dktl, 2805 "dktregions":dktp, 2806 "dktcortex":dktc, 2807 "bf":deep_bf['description'], 2808 } 2809 2810 outputs = { 2811 "brain_n4_dnz": img, 2812 "brain_extraction": imgbxt, 2813 "left_right": mylr, 2814 "dkt_parc": myparc, 2815 "bf":deep_bf['segmentation'], 2816 "dataframes": mydataframes 2817 } 2818 2819 return outputs 2820 2821 2822 ##### traditional deformable registration as a high-resolution complement to above 2823 wm_tractsL = None 2824 wm_tractsR = None 2825 wmtdfL = None 2826 wmtdfR = None 2827 reg = None 2828 if labels_to_register is not None and not is_test: 2829 reg = hemi_reg( 2830 input_image = img, 2831 input_image_tissue_segmentation = myparc['tissue_segmentation'], 2832 input_image_hemisphere_segmentation = mylr, 2833 input_template=templatea, 2834 input_template_hemisphere_labels=templatealr, 2835 output_prefix = output_prefix + "_SYN", 2836 labels_to_register = labels_to_register, 2837 is_test=is_test ) 2838 if verbose: 2839 print("wm tracts") 2840 ##### how to use the hemi-reg output to generate any roi value from a template roi 2841 wm_tracts = ants.image_read( get_data( "wm_major_tracts", target_extension='.nii.gz' ) ) 2842 wm_tractsL = ants.apply_transforms( img, wm_tracts, reg['synL']['invtransforms'], 2843 interpolator='nearestNeighbor' ) * ants.threshold_image( mylr, 1, 1 ) 2844 wm_tractsR = ants.apply_transforms( img, wm_tracts, reg['synR']['invtransforms'], 2845 interpolator='nearestNeighbor' ) * ants.threshold_image( mylr, 2, 2 ) 2846 wmtdfL = map_segmentation_to_dataframe( "wm_major_tracts", wm_tractsL ) 2847 wmtdfR = map_segmentation_to_dataframe( "wm_major_tracts", wm_tractsR ) 2848 2849 if verbose: 2850 print("hippocampus") 2851 2852 ##### specialized labeling for hippocampus 2853 do_hippmapper = False 2854 if do_hippmapper: 2855 ntries = 10 2856 if is_test: 2857 ntries = 1 2858 hippLR = deep_hippo( img=img, template=templateb, number_of_tries=ntries, 2859 tx_type='Affine', sr_model = sr_model ) 2860 2861 if verbose: 2862 print("medial temporal lobe") 2863 2864 ##### deep_flash medial temporal lobe parcellation 2865 deep_flash = deep_mtl(img, sr_model = sr_model ) 2866 2867 if verbose: 2868 print("deep CIT168") 2869 ##### deep CIT168 segmentation - relatively fast 2870 deep_cit = deep_cit168( img, priors = cit168lab, 2871 binary_mask = braintissuemask ) 2872 2873 if verbose: 2874 print( "SN-specific segmentation" ) 2875# input_image_region_segmentation, input_template, input_template_region_segmentation, output_prefix, padding=10, labels_to_register=[2, 3, 4, 5], total_sigma=0.5, is_test=False) 2876 2877 if not is_test: 2878 tbinseg = ants.mask_image( cit168labT, cit168labT, [7,9,23,25,33,34], binarize=True) 2879 tbinseg = ants.morphology( tbinseg, "dilate", 14 ) 2880 ibinseg = ants.apply_transforms( img, tbinseg, cit168reg['invtransforms'], 2881 interpolator='nearestNeighbor') 2882 snreg = region_reg( img, myparc['tissue_segmentation'], ibinseg, 2883 cit168adni, tbinseg, output_prefix=output_prefix + "_SNREG", 2884 padding = 4, is_test=False )['synL'] 2885 tbinseg = ants.mask_image( cit168labT, cit168labT, [7,9,23,25,33,34], binarize=False) 2886 snseg = ants.apply_transforms( img, tbinseg, 2887 snreg['invtransforms'], interpolator = 'nearestNeighbor' ) 2888 snseg = snseg * ants.threshold_image( myparc['tissue_segmentation'], 2, 6 ) 2889 snseg_desc = map_segmentation_to_dataframe( 'CIT168_Reinf_Learn_v1_label_descriptions_pad', snseg ).dropna(axis=0) 2890 else : 2891 snseg = None 2892 snseg_desc = None 2893 2894 if verbose: 2895 print( "brainstem and cerebellar segmentation" ) 2896 # midbrain/brainstem 2897 brainstemseg = ants.apply_transforms( img, cit168labStem, 2898 cit168reg['invtransforms'], interpolator = 'nearestNeighbor' ) 2899 brainstemseg = brainstemseg * braintissuemask 2900 brainstem_desc = map_segmentation_to_dataframe( 'CIT168_T1w_700um_pad_adni_brainstem', brainstemseg ) 2901 brainstem_desc = brainstem_desc.loc[:, ~brainstem_desc.columns.str.contains('^Side')] 2902 2903 2904 # cerebellum 2905 cereb = antspynet.cerebellum_morphology( ants.iMath( x, "Normalize" ), compute_thickness_image=False, verbose=False, do_preprocessing=True ) 2906 # refinement with cerebellum estimate (comment out since it's not needed for this image). 2907 maskc = ants.threshold_image(cereb['cerebellum_probability_image'], 0.5, 1, 1, 0) 2908 cereb = antspynet.cerebellum_morphology( ants.iMath( x, "Normalize" ), cerebellum_mask=maskc, compute_thickness_image=False, verbose=False, do_preprocessing=True ) 2909 cereb_desc = map_segmentation_to_dataframe( 'cerebellum', cereb['parcellation_segmentation_image'] ).dropna(axis=0) 2910 2911 if verbose: 2912 print( "antspyt1w.hierarchical complete" ) 2913 mydataframes = { 2914 "rbp": myqc['brain'], 2915 "icv": myicv, 2916 "hemispheres":hemi, 2917 "tissues":tissue, 2918 "dktlobes":dktl, 2919 "dktregions":dktp, 2920 "dktcortex":dktc, 2921 "wmtracts_left":wmtdfL, 2922 "wmtracts_right":wmtdfR, 2923# "wmh":myhypo['wmh_summary'], 2924 "mtl":deep_flash['mtl_description'], 2925 "bf":deep_bf['description'], 2926 "cit168":cit168lab_desc, 2927 "deep_cit168":deep_cit['description'], 2928 "snseg":snseg_desc, 2929# "hippLR":hippLR['description'], 2930 "brainstem": brainstem_desc, 2931 "cerebellum": cereb_desc 2932 } 2933 2934 outputs = { 2935 "brain_n4_dnz": img, 2936 "brain_n4_dnz_png": myqc['brain_image'], 2937 "brain_extraction": imgbxt, 2938 "tissue_seg_png": tissue_seg_png, 2939 "left_right": mylr, 2940 "dkt_parc": myparc, 2941 "registration":reg, 2942# "hippLR":hippLR['segmentation'], 2943# "white_matter_hypointensity":myhypo['wmh_probability_image'], 2944 "wm_tractsL":wm_tractsL, 2945 "wm_tractsR":wm_tractsR, 2946 "mtl":deep_flash['mtl_segmentation'], 2947 "bf":deep_bf['segmentation'], 2948 "deep_cit168lab": deep_cit['segmentation'], 2949 "cit168lab": cit168lab, 2950 "cit168reg": cit168reg, 2951 "snseg":snseg, 2952 "snreg":snreg, 2953 "brainstem": brainstemseg, 2954 "cerebellum":cereb['parcellation_segmentation_image'], 2955 "dataframes": mydataframes 2956 } 2957 2958 return outputs
Default processing for a T1-weighted image. See README.
Arguments
x : T1-weighted neuroimage antsImage
output_prefix: string directory and prefix
labels_to_register: list of integer segmentation labels (of 1 to 6 as defined by atropos: csf, gm, wm, dgm, brainstem, cerebellum) to define the tissue types / regions of the brain to register. set to None to skip registration which will be faster but omit some results.
imgbxt : pre-existing brain extraction - a binary image - will disable some processing
img6seg: pre-existing brain segmentation - a 6-class image ( ANTs standard )
cit168 : boolean returns labels from CIT168 atlas with high-resolution registration otherwise, low-resolution regitration is used.
is_test: boolean ( parameters to run more quickly but with low quality )
atropos_prior: prior weight for atropos post-processing; set to None if you do not want to use this. will modify the CSF, GM and WM segmentation to better fit the image intensity at the resolution of the input image.
sr_model: optional will trigger sr-based upsampling in some functions.
verbose: boolean
Returns
dataframes and associated derived data
- brain_n4_dnz : extracted brain denoised and bias corrected
- brain_extraction : brain mask
- rbp: random basis projection results
- left_right : left righ hemisphere segmentation
- dkt_parc : dictionary object containing segmentation labels
- registration : dictionary object containing registration results
- hippLR : dictionary object containing hippocampus results
- medial_temporal_lobe : dictionary object containing deep_flash (medial temporal lobe parcellation) results
- white_matter_hypointensity : dictionary object containing WMH results
- wm_tractsL : white matter tracts, left
- wm_tractsR : white matter tracts, right
- dataframes : summary data frames
def
random_basis_projection( x, template, type_of_transform='Similarity', refbases=None, nBasis=10, random_state=99):
477def random_basis_projection( x, template, 478 type_of_transform='Similarity', 479 refbases = None, 480 nBasis=10, random_state = 99 ): 481 """ 482 Produce unbiased data descriptors for a given image which can be used 483 to assist data inspection and ranking. can be used with any image 484 brain extracted or not, any modality etc. but we assume we can 485 meaningfully map to a template, at least with a low-dimensional 486 transformation, e.g. Translation, Rigid, Similarity. 487 488 Arguments 489 --------- 490 x : antsImage 491 492 template : antsImage reference template 493 494 type_of_transform: one of Translation, Rigid, Similarity, Affine 495 496 refbases : reference bases for outlierness calculations 497 498 nBasis : number of variables to derive 499 500 random_state : seed 501 502 Returns 503 ------- 504 dataframe with projections and an outlierness estimate. 505 506 the outlierness estimate is based on a small reference dataset of young controls. 507 the user/researcher may want to use a different reference set. see the 508 function loop_outlierness for one way to do that. 509 510 """ 511 template = ants.crop_image( template ) 512 template = ants.iMath( template, "Normalize" ) 513 np.random.seed(int(random_state)) 514 nvox = template.shape 515 # X = np.random.rand( nBasis+1, myproduct( nvox ) ) 516 # u, s, randbasis = svds(X, k=nBasis) 517 # if randbasis.shape[1] != myproduct(nvox): 518 # raise ValueError("columns in rand basis do not match the nvox product") 519 520 randbasis = np.random.randn( myproduct( nvox ), nBasis ) 521 rbpos = randbasis.copy() 522 rbpos[rbpos<0] = 0 523 norm = ants.iMath( x, "Normalize" ) 524 trans = ants.registration( template, norm, 525 type_of_transform='antsRegistrationSyNQuickRepro[t]' ) 526 resamp = ants.registration( template, norm, 527 type_of_transform=type_of_transform, total_sigma=0.5, 528 # aff_metric='GC', 529 random_seed=1, initial_transform=trans['fwdtransforms'][0] )['warpedmovout'] 530# mydelta = ants.from_numpy( ( resamp - template ).abs() ) 531 mydelta = resamp - template 532 imat = ants.get_neighborhood_in_mask( mydelta, mydelta*0+1,[0,0,0], boundary_condition='mean' ) 533 uproj = np.matmul(imat, randbasis) 534 uprojpos = np.matmul(imat, rbpos) 535 record = {} 536 uproj_counter = 0 537 for i in uproj[0]: 538 uproj_counter += 1 539 name = "RandBasisProj" + str(uproj_counter).zfill(2) 540 record[name] = i 541 uprojpos_counter = 0 542 for i in uprojpos[0]: 543 uprojpos_counter += 1 544 name = "RandBasisProjPos" + str(uprojpos_counter).zfill(2) 545 record[name] = i 546 df = pd.DataFrame(record, index=[0]) 547 548 if refbases is None: 549 refbases = pd.read_csv( get_data( "reference_basis", target_extension='.csv' ) ) 550 df['loop_outlier_probability'] = loop_outlierness( df, refbases, 551 n_neighbors=refbases.shape[0] )[ refbases.shape[0] ] 552 mhdist = 0.0 553 if nBasis == 10: 554 temp = pd.concat( [ refbases, df.iloc[:,:nBasis] ], axis=0 ) 555 mhdist = mahalanobis_distance( temp )['distance'][ refbases.shape[0] ] 556 df['mhdist'] = mhdist 557 df['templateL1']=mydelta.abs().mean() 558 return df
Produce unbiased data descriptors for a given image which can be used to assist data inspection and ranking. can be used with any image brain extracted or not, any modality etc. but we assume we can meaningfully map to a template, at least with a low-dimensional transformation, e.g. Translation, Rigid, Similarity.
Arguments
x : antsImage
template : antsImage reference template
type_of_transform: one of Translation, Rigid, Similarity, Affine
refbases : reference bases for outlierness calculations
nBasis : number of variables to derive
random_state : seed
Returns
dataframe with projections and an outlierness estimate.
the outlierness estimate is based on a small reference dataset of young controls. the user/researcher may want to use a different reference set. see the function loop_outlierness for one way to do that.
deep_dkt
def
deep_hippo( img, template, number_of_tries=10, tx_type='antsRegistrationSyNQuickRepro[a]', sr_model=None, verbose=False):
1183def deep_hippo( 1184 img, 1185 template, 1186 number_of_tries = 10, 1187 tx_type="antsRegistrationSyNQuickRepro[a]", 1188 sr_model = None, 1189 verbose=False 1190): 1191 1192 1193 super_resolution = None 1194 1195 avgleft = img * 0 1196 avgright = img * 0 1197 for k in range(number_of_tries): 1198 if verbose: 1199 print( "try " + str(k)) 1200 rig = ants.registration( template, ants.rank_intensity(img), 1201 tx_type, random_seed=k, total_sigma=0.5, verbose=verbose ) 1202 if verbose: 1203 print( rig ) 1204 rigi = ants.apply_transforms( template, img, rig['fwdtransforms'] ) 1205 if verbose: 1206 print( "done with warp - do hippmapp3r" ) 1207 hipp = antspynet.hippmapp3r_segmentation( rigi, do_preprocessing=False ) 1208 if verbose: 1209 print( "done with hippmapp3r - backtransform" ) 1210 1211 if sr_model is not None: 1212 mysr = super_resolution_segmentation_per_label( 1213 rigi, segmentation=hipp, upFactor=[2,2,2], 1214 sr_model=sr_model, segmentation_numbers=[1,2], 1215 dilation_amount=0, 1216 probability_images=None, 1217 probability_labels=[1,2], max_lab_plus_one=True, verbose=True 1218 ) 1219 hipp = mysr['super_resolution_segmentation'] 1220 1221 hippr = ants.apply_transforms( 1222 img, 1223 hipp, 1224 rig['fwdtransforms'], 1225 whichtoinvert=[True], 1226 interpolator='nearestNeighbor', 1227 ) 1228 avgleft = avgleft + ants.threshold_image( hippr, 2, 2 ).iMath("GetLargestComponent",1) / float(number_of_tries) 1229 avgright = avgright + ants.threshold_image( hippr, 1, 1 ).iMath("GetLargestComponent",1) / float(number_of_tries) 1230 1231 1232 avgright = ants.iMath(avgright,"Normalize") # output: probability image right 1233 avgleft = ants.iMath(avgleft,"Normalize") # output: probability image left 1234 hippright_bin = ants.threshold_image( avgright, 0.5, 2.0 ).iMath("GetLargestComponent",1) 1235 hippleft_bin = ants.threshold_image( avgleft, 0.5, 2.0 ).iMath("GetLargestComponent",1) 1236 hipp_bin = hippleft_bin + hippright_bin * 2 1237 1238 hippleftORlabels = ants.label_geometry_measures(hippleft_bin, avgleft) 1239 hippleftORlabels['Description'] = 'left hippocampus' 1240 hipprightORlabels = ants.label_geometry_measures(hippright_bin, avgright) 1241 hipprightORlabels['Description'] = 'right hippocampus' 1242 # hippleftORlabels=hippleftORlabels.append( hipprightORlabels ) 1243 # FIXME possible fix below 1244 hippleftORlabels = pd.concat( [hippleftORlabels, hipprightORlabels] , axis=1 ) 1245 hippleftORlabels['Label']=[1,2] 1246 labels = { 1247 'segmentation':hipp_bin, 1248 'description':hippleftORlabels, 1249 'HLProb':avgleft, 1250 'HRProb':avgright 1251 } 1252 return labels
def
deep_tissue_segmentation( x, template=None, registration_map=None, atropos_prior=None, sr_model=None):
1031def deep_tissue_segmentation( x, template=None, registration_map=None, 1032 atropos_prior=None, sr_model=None ): 1033 """ 1034 modified slightly more efficient deep atropos that also handles the 1035 extra CSF issue. returns segmentation and probability images. see 1036 the tissues csv available from get_data. 1037 1038 x: input image, raw 1039 1040 template: MNI space template, should be "croppedMni152" or "biobank" 1041 1042 registration_map: pre-existing output from ants.registration 1043 1044 atropos_prior: prior weight for atropos post-processing; set to None if you 1045 do not want to use this. will modify the CSF, GM and WM segmentation to 1046 better fit the image intensity at the resolution of the input image. 1047 1048 sr_model: optional (FIXME) 1049 1050 """ 1051 1052 dapper = antspynet.deep_atropos( x, do_denoising=False ) 1053 1054 myk='segmentation_image' 1055 if atropos_prior is not None: 1056 msk = ants.threshold_image( dapper['segmentation_image'], 2, 3 ).iMath("GetLargestComponent",50) 1057 msk = ants.morphology( msk, "close", 2 ) 1058 mskfull = ants.threshold_image( dapper['segmentation_image'], 1, 6 ) 1059 mskfull = mskfull - msk 1060 priors = dapper[myk][1:4] 1061 for k in range( len( priors ) ): 1062 priors[k]=ants.image_clone( priors[k]*msk ) 1063 aap = ants.atropos( x, msk, i=priors, m='[0.0,1x1x1]', 1064 c = '[1,0]', priorweight=atropos_prior, verbose=1 ) 1065 dapper['segmentation_image'] = aap['segmentation'] * msk + dapper['segmentation_image'] * mskfull 1066 dapper['probability_images'][1:4] = aap['probabilityimages'] 1067 1068 return dapper
modified slightly more efficient deep atropos that also handles the extra CSF issue. returns segmentation and probability images. see the tissues csv available from get_data.
x: input image, raw
template: MNI space template, should be "croppedMni152" or "biobank"
registration_map: pre-existing output from ants.registration
atropos_prior: prior weight for atropos post-processing; set to None if you do not want to use this. will modify the CSF, GM and WM segmentation to better fit the image intensity at the resolution of the input image.
sr_model: optional (FIXME)
def
deep_brain_parcellation( target_image, template, img6seg=None, do_cortical_propagation=False, atropos_prior=None, verbose=True):
1070def deep_brain_parcellation( 1071 target_image, 1072 template, 1073 img6seg = None, 1074 do_cortical_propagation=False, 1075 atropos_prior=None, 1076 verbose=True, 1077): 1078 """ 1079 modified slightly more efficient deep dkt that also returns atropos output 1080 thus providing a complete hierarchical parcellation of t1w. we run atropos 1081 here so we dont need to redo registration separately. see 1082 the lobes and dkt csv available from get_data. 1083 1084 target_image: input image 1085 1086 template: MNI space template, should be "croppedMni152" or "biobank" 1087 1088 img6seg: optional pre-existing brain segmentation - a 6-class image ( ANTs standard ) 1089 1090 do_cortical_propagation: boolean, adds a bit extra time to propagate cortical 1091 labels explicitly into cortical segmentation 1092 1093 atropos_prior: prior weight for atropos post-processing; set to None if you 1094 do not want to use this. will modify the CSF, GM and WM segmentation to 1095 better fit the image intensity at the resolution of the input image. 1096 1097 verbose: boolean 1098 1099 1100 Returns 1101 ------- 1102 a dictionary containing: 1103 1104 - tissue_segmentation : 6 tissue segmentation 1105 - tissue_probabilities : probability images associated with above 1106 - dkt_parcellation : tissue agnostic DKT parcellation 1107 - dkt_lobes : major lobes of the brain 1108 - dkt_cortex: cortical tissue DKT parcellation (if requested) 1109 - hemisphere_labels: free to get hemisphere labels 1110 - wmSNR : white matter signal-to-noise ratio 1111 - wmcsfSNR : white matter to csf signal-to-noise ratio 1112 1113 """ 1114 if verbose: 1115 print("Begin registration") 1116 1117 rig = ants.registration( template, ants.rank_intensity(target_image), 1118 "antsRegistrationSyNQuickRepro[a]", 1119 aff_iterations = (500,200,0,0), 1120 total_sigma=0.5, 1121 random_seed=1 ) 1122 rigi = ants.apply_transforms( template, target_image, rig['fwdtransforms']) 1123 1124 if verbose: 1125 print("Begin Atropos tissue segmentation") 1126 1127 if img6seg is None: 1128 if verbose: 1129 print( "do deep_tissue_segmentation") 1130 mydap = deep_tissue_segmentation( target_image ) 1131 else: 1132 if verbose: 1133 print( "use existing deep_tissue_segmentation") 1134 mydap = { 'segmentation_image': img6seg, 'probability_images': None } 1135 1136 if verbose: 1137 print("End Atropos tissue segmentation") 1138 1139 if verbose: 1140 print("Begin DKT") 1141 1142 dktprepro = True 1143 dkt = antspynet.desikan_killiany_tourville_labeling( 1144 target_image, 1145 do_preprocessing=True, 1146 return_probability_images=False, 1147 do_lobar_parcellation = True, 1148 do_denoising=False 1149 ) 1150 1151 myhemiL = ants.threshold_image( dkt['lobar_parcellation'], 1, 6 ) 1152 myhemiR = ants.threshold_image( dkt['lobar_parcellation'], 7, 12 ) 1153 myhemi = myhemiL + myhemiR * 2.0 1154 brainmask = ants.threshold_image( mydap['segmentation_image'], 1, 6 ) 1155 myhemi = ants.iMath( brainmask, 'PropagateLabelsThroughMask', myhemi, 100, 0) 1156 1157 cortprop = None 1158 if do_cortical_propagation: 1159 cortprop = ants.threshold_image( mydap['segmentation_image'], 2, 2 ) 1160 cortlab = dkt['segmentation_image'] * ants.threshold_image( dkt['segmentation_image'], 1000, 5000 ) 1161 cortprop = ants.iMath( cortprop, 'PropagateLabelsThroughMask', 1162 cortlab, 1, 0) 1163 1164 wmseg = ants.threshold_image( mydap['segmentation_image'], 3, 3 ) 1165 wmMean = target_image[ wmseg == 1 ].mean() 1166 wmStd = target_image[ wmseg == 1 ].std() 1167 csfseg = ants.threshold_image( mydap['segmentation_image'], 1, 1 ) 1168 csfStd = target_image[ csfseg == 1 ].std() 1169 wmSNR = wmMean/wmStd 1170 wmcsfSNR = wmMean/csfStd 1171 1172 return { 1173 "tissue_segmentation":mydap['segmentation_image'], 1174 "tissue_probabilities":mydap['probability_images'], 1175 "dkt_parcellation":dkt['segmentation_image'], 1176 "dkt_lobes":dkt['lobar_parcellation'], 1177 "dkt_cortex": cortprop, 1178 "hemisphere_labels": myhemi, 1179 "wmSNR": wmSNR, 1180 "wmcsfSNR": wmcsfSNR, }
modified slightly more efficient deep dkt that also returns atropos output thus providing a complete hierarchical parcellation of t1w. we run atropos here so we dont need to redo registration separately. see the lobes and dkt csv available from get_data.
target_image: input image
template: MNI space template, should be "croppedMni152" or "biobank"
img6seg: optional pre-existing brain segmentation - a 6-class image ( ANTs standard )
do_cortical_propagation: boolean, adds a bit extra time to propagate cortical labels explicitly into cortical segmentation
atropos_prior: prior weight for atropos post-processing; set to None if you do not want to use this. will modify the CSF, GM and WM segmentation to better fit the image intensity at the resolution of the input image.
verbose: boolean
Returns
a dictionary containing:
- tissue_segmentation : 6 tissue segmentation
- tissue_probabilities : probability images associated with above
- dkt_parcellation : tissue agnostic DKT parcellation
- dkt_lobes : major lobes of the brain
- dkt_cortex: cortical tissue DKT parcellation (if requested)
- hemisphere_labels: free to get hemisphere labels
- wmSNR : white matter signal-to-noise ratio
- wmcsfSNR : white matter to csf signal-to-noise ratio
def
deep_mtl(t1, sr_model=None, verbose=True):
1384def deep_mtl(t1, sr_model=None, verbose=True): 1385 1386 """ 1387 Hippocampal/Enthorhinal segmentation using "Deep Flash" 1388 1389 Perform hippocampal/entorhinal segmentation in T1 images using 1390 labels from Mike Yassa's lab 1391 1392 https://faculty.sites.uci.edu/myassa/ 1393 1394 The labeling is as follows: 1395 Label 0 : background 1396 Label 5 : left aLEC 1397 Label 6 : right aLEC 1398 Label 7 : left pMEC 1399 Label 8 : right pMEC 1400 Label 9 : left perirhinal 1401 Label 10: right perirhinal 1402 Label 11: left parahippocampal 1403 Label 12: right parahippocampal 1404 Label 13: left DG/CA3 1405 Label 14: right DG/CA3 1406 Label 15: left CA1 1407 Label 16: right CA1 1408 Label 17: left subiculum 1409 Label 18: right subiculum 1410 1411 """ 1412 1413 verbose = False 1414 1415 labels = (0, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18) 1416 label_descriptions = ['background', 1417 'left aLEC', 1418 'right aLEC', 1419 'left pMEC', 1420 'right pMEC', 1421 'left perirhinal', 1422 'right perirhinal', 1423 'left parahippocampal', 1424 'right parahippocampal', 1425 'left DG/CA3', 1426 'right DG/CA3', 1427 'left CA1', 1428 'right CA1', 1429 'left subiculum', 1430 'right subiculum' 1431 ] 1432 1433 template = ants.image_read(antspynet.get_antsxnet_data("deepFlashTemplateT1SkullStripped")) 1434 registration = ants.registration(fixed=template, moving=t1, 1435 type_of_transform="antsRegistrationSyNQuickRepro[a]", total_sigma=0.5, verbose=verbose) 1436 template_transforms = dict(fwdtransforms=registration['fwdtransforms'], 1437 invtransforms=registration['invtransforms']) 1438 t1_warped = registration['warpedmovout'] 1439 1440 df = antspynet.deep_flash(t1_warped, 1441 do_preprocessing=False, 1442 use_rank_intensity = True, 1443 verbose=verbose) 1444 1445 if False: 1446 if sr_model is not None: 1447 ulabs = np.unique( df['segmentation_image'].numpy() ) 1448 ulabs.sort() 1449 ulabs = ulabs[1:len(ulabs)] 1450 ulabs = ulabs.tolist() 1451 newprobs = super_resolution_segmentation_per_label( 1452 t1_warped, 1453 df['segmentation_image'], 1454 [2,2,2], 1455 sr_model, 1456 ulabs, 1457 probability_images=df['probability_images'][1:len(df['probability_images'])], 1458 probability_labels=ulabs, 1459 max_lab_plus_one=True ) 1460 for k in range(1, len(df['probability_images']) ): 1461 df['probability_images'][k] = newprobs['probability_images'][k-1] 1462 1463 if sr_model is not None: 1464 newprobs = super_resolution_segmentation_with_probabilities( 1465 t1_warped, df['probability_images'], sr_model ) 1466 df['probability_images'] = newprobs['sr_probabilities'] 1467 1468 if False: 1469 probability_images = list() 1470 for i in range(len(df['probability_images'])): 1471 probability_image = ants.apply_transforms(fixed=t1, 1472 moving=df['probability_images'][i], 1473 transformlist=template_transforms['invtransforms'], 1474 whichtoinvert=[True], 1475 interpolator="linear", 1476 verbose=verbose) 1477 probability_images.append(probability_image) 1478 image_matrix = ants.image_list_to_matrix(probability_images[1:(len(probability_images))], t1 * 0 + 1) 1479 background_foreground_matrix = np.stack([ants.image_list_to_matrix([probability_images[0]], t1 * 0 + 1), 1480 np.expand_dims(np.sum(image_matrix, axis=0), axis=0)]) 1481 foreground_matrix = np.argmax(background_foreground_matrix, axis=0) 1482 segmentation_matrix = (np.argmax(image_matrix, axis=0) + 1) * foreground_matrix 1483 segmentation_image = ants.matrix_to_images(np.expand_dims(segmentation_matrix, axis=0), t1 * 0 + 1)[0] 1484 relabeled_image = ants.image_clone(segmentation_image) 1485 for i in range(len(labels)): 1486 relabeled_image[segmentation_image==i] = labels[i] 1487 1488 relabeled_image = ants.apply_transforms( fixed=t1, 1489 moving=df['segmentation_image'], 1490 transformlist=template_transforms['invtransforms'], 1491 whichtoinvert=[True], 1492 interpolator="nearestNeighbor", 1493 verbose=verbose) 1494 mtl_description = map_segmentation_to_dataframe( 'mtl_description', relabeled_image ) 1495 1496 deep_mtl_dictionary = { 1497 'mtl_description':mtl_description, 1498 'mtl_segmentation':relabeled_image 1499# 'mtl_probability_images':probability_images 1500 } 1501 return(deep_mtl_dictionary)
Hippocampal/Enthorhinal segmentation using "Deep Flash"
Perform hippocampal/entorhinal segmentation in T1 images using labels from Mike Yassa's lab
https://faculty.sites.uci.edu/myassa/
The labeling is as follows: Label 0 : background Label 5 : left aLEC Label 6 : right aLEC Label 7 : left pMEC Label 8 : right pMEC Label 9 : left perirhinal Label 10: right perirhinal Label 11: left parahippocampal Label 12: right parahippocampal Label 13: left DG/CA3 Label 14: right DG/CA3 Label 15: left CA1 Label 16: right CA1 Label 17: left subiculum Label 18: right subiculum
def
label_hemispheres(x, template=None, templateLR=None, reg_iterations=[200, 50, 2, 0]):
994def label_hemispheres( x, template=None, templateLR=None, reg_iterations=[200,50,2,0] ): 995 """ 996 quick somewhat noisy registration solution to hemisphere labeling. typically 997 we label left as 1 and right as 2. 998 999 x: input image 1000 1001 template: MNI space template, should be "croppedMni152" or "biobank" 1002 1003 templateLR: a segmentation image of template hemispheres 1004 1005 reg_iterations: reg_iterations for ants.registration 1006 1007 """ 1008 1009 if template is None or templateLR is None: 1010 tfn = get_data('T_template0', target_extension='.nii.gz' ) 1011 tfnw = get_data('T_template0_WMP', target_extension='.nii.gz' ) 1012 tlrfn = get_data('T_template0_LR', target_extension='.nii.gz' ) 1013 bfn = antspynet.get_antsxnet_data( "croppedMni152" ) 1014 template = ants.image_read( tfn ) 1015 template = ( template * antspynet.brain_extraction( template, 't1' ) ).iMath( "Normalize" ) 1016 templateawmprior = ants.image_read( tfnw ) 1017 templateLR = ants.image_read( tlrfn ) 1018 1019 reg = ants.registration( 1020 ants.rank_intensity(x), 1021 ants.rank_intensity(template), 1022 'antsRegistrationSyNQuickRepro[s]', 1023 aff_metric='GC', 1024 syn_metric='CC', 1025 syn_sampling=2, 1026 reg_iterations=reg_iterations, total_sigma=0.5, 1027 random_seed = 1 ) 1028 return( ants.apply_transforms( x, templateLR, reg['fwdtransforms'], 1029 interpolator='nearestNeighbor') )
quick somewhat noisy registration solution to hemisphere labeling. typically we label left as 1 and right as 2.
x: input image
template: MNI space template, should be "croppedMni152" or "biobank"
templateLR: a segmentation image of template hemispheres
reg_iterations: reg_iterations for ants.registration
def
brain_extraction(x, dilation=8.0, method='v1', deform=True, verbose=False):
794def brain_extraction( x, dilation = 8.0, method = 'v1', deform=True, verbose=False ): 795 """ 796 quick brain extraction for individual images 797 798 x: input image 799 800 dilation: amount to dilate first brain extraction in millimeters 801 802 method: version currently v0 or any other string gives two different results 803 804 deform: map the image to the training hypersphere before running the extraction 805 806 verbose: boolean 807 808 """ 809 return antspynet.brain_extraction( ants.iMath( x, "Normalize" ), modality="t1threetissue")['segmentation_image'].threshold_image(1,1) 810 811 if deform: 812 reorient_template_file_name_path = antspynet. get_antsxnet_data("S_template3" ) 813 template = ants.image_read( reorient_template_file_name_path ) 814 reg = ants.registration( template, x, 'antsRegistrationSyNQuickRepro[s]', 815 total_sigma=0.5 ) 816 817 closedilmm = 5.0 818 spacing = ants.get_spacing(x) 819 spacing_product = spacing[0] * spacing[1] * spacing[2] 820 spcmin = min( spacing ) 821 dilationRound = int(np.round( dilation / spcmin )) 822 closedilRound = int(np.round( closedilmm / spcmin )) 823 if deform: 824 xn3 = ants.n3_bias_field_correction( reg['warpedmovout'], 8 ).n3_bias_field_correction( 4 ) 825 xn3 = ants.iMath(xn3, "TruncateIntensity",0.001,0.999).iMath("Normalize") 826 else: 827 xn3 = ants.n3_bias_field_correction( x, 8 ).n3_bias_field_correction( 4 ) 828 xn3 = ants.iMath(xn3, "TruncateIntensity",0.001,0.999).iMath("Normalize") 829 if method == 'v0': 830 if verbose: 831 print("method v0") 832 bxtmethod = 't1combined[' + str(closedilRound) +']' # better for individual subjects 833 bxt = antspynet.brain_extraction( xn3, bxtmethod ).threshold_image(2,3).iMath("GetLargestComponent",50).iMath("FillHoles") 834 if deform: 835 bxt = ants.apply_transforms( x, bxt, reg['invtransforms'], interpolator='nearestNeighbor' ) 836 return bxt 837 if method == 'v1': 838 if verbose: 839 print("method v1") 840 bxt = antspynet.brain_extraction( xn3, 't1' ).threshold_image(0.5,3).iMath("GetLargestComponent",50).iMath("FillHoles") 841 if deform: 842 bxt = ants.apply_transforms( x, bxt, reg['invtransforms'], interpolator='nearestNeighbor' ) 843 return bxt 844 if verbose: 845 print("method candidate") 846 bxt0 = antspynet.brain_extraction( xn3, "t1" ).threshold_image(0.5,1.0).iMath("GetLargestComponent",50).morphology( "close", closedilRound ).iMath("FillHoles") 847 bxt0dil = ants.iMath( bxt0, "MD", dilationRound ) 848 image = ants.iMath( xn3 * bxt0dil,"Normalize")*255 849 # no no brainer 850 model = antspynet.create_nobrainer_unet_model_3d((None, None, None, 1)) 851 weights_file_name = antspynet.get_pretrained_network("brainExtractionNoBrainer") 852 model.load_weights(weights_file_name) 853 image_array = image.numpy() 854 image_robust_range = np.quantile(image_array[np.where(image_array != 0)], (0.02, 0.98)) 855 threshold_value = 0.10 * (image_robust_range[1] - image_robust_range[0]) + image_robust_range[0] 856 thresholded_mask = ants.threshold_image(image, -10000, threshold_value, 0, 1) 857 thresholded_image = image * thresholded_mask 858 image_resampled = ants.resample_image(thresholded_image, (256, 256, 256), use_voxels=True) 859 image_array = np.expand_dims(image_resampled.numpy(), axis=0) 860 image_array = np.expand_dims(image_array, axis=-1) 861 brain_mask_array = np.squeeze(model.predict(image_array)) 862 brain_mask_resampled = ants.copy_image_info(image_resampled, ants.from_numpy(brain_mask_array)) 863 brain_mask_image = ants.resample_image(brain_mask_resampled, image.shape, use_voxels=True, interp_type=1) 864 brain_mask_image = brain_mask_image * ants.threshold_image( brain_mask_image, 0.5, 1e9 ) 865 minimum_brain_volume = round(649933.7/spacing_product) 866 bxt1 = ants.label_clusters(brain_mask_image, minimum_brain_volume) 867 nblabels = np.unique( bxt1.numpy() )# .astype(int) 868 maxol = 0.0 869 bestlab = bxt0 870 871 for j in range(1,len(nblabels)): 872 temp = ants.threshold_image( bxt1, j, j ) 873 tempsum = ( bxt0 * temp ).sum() 874 dicer = ants.label_overlap_measures( temp, bxt0 ) 875 if tempsum > maxol and dicer['MeanOverlap'][1] > 0.5 : 876 if verbose: 877 print( str(j) + ' dice ' + str( dicer['MeanOverlap'][1] ) ) 878 maxol = tempsum 879 bestlab = temp 880 adder = trim_segmentation_by_distance( bxt0, 1, closedilmm ) 881 bestlab = ants.threshold_image( bestlab + adder, 1, 2 ) 882 return bestlab
quick brain extraction for individual images
x: input image
dilation: amount to dilate first brain extraction in millimeters
method: version currently v0 or any other string gives two different results
deform: map the image to the training hypersphere before running the extraction
verbose: boolean
def
hemi_reg( input_image, input_image_tissue_segmentation, input_image_hemisphere_segmentation, input_template, input_template_hemisphere_labels, output_prefix, padding=10, labels_to_register=[2, 3, 4, 5], total_sigma=0.5, is_test=False):
1607def hemi_reg( 1608 input_image, 1609 input_image_tissue_segmentation, 1610 input_image_hemisphere_segmentation, 1611 input_template, 1612 input_template_hemisphere_labels, 1613 output_prefix, 1614 padding=10, 1615 labels_to_register = [2,3,4,5], 1616 total_sigma=0.5, 1617 is_test=False ): 1618 """ 1619 hemisphere focused registration that will produce jacobians and figures to 1620 support data inspection 1621 1622 input_image: input image 1623 1624 input_image_tissue_segmentation: segmentation produced in ANTs style ie with 1625 labels defined by atropos brain segmentation (1 to 6) 1626 1627 input_image_hemisphere_segmentation: left (1) and right (2) hemisphere 1628 segmentation 1629 1630 input_template: template to which we register; prefer a population-specific 1631 relatively high-resolution template instead of MNI or biobank. 1632 1633 input_template_hemisphere_labels: a segmentation image of template hemispheres 1634 with left labeled 1 and right labeled 2 1635 1636 output_prefix: a path and prefix for registration related outputs 1637 1638 padding: number of voxels to pad images, needed for diffzero 1639 1640 labels_to_register: list of integer segmentation labels to use to define 1641 the tissue types / regions of the brain to register. 1642 1643 total_sigma: scalar >= 0.0 ; higher means more constrained registration. 1644 1645 is_test: boolean. this function can be long running by default. this would 1646 help testing more quickly by running fewer iterations. 1647 1648 """ 1649 1650 img = ants.rank_intensity( input_image ) 1651 # ionlycerebrum = ants.mask_image( input_image_tissue_segmentation, 1652 # input_image_tissue_segmentation, labels_to_register, 1 ) 1653 ionlycerebrum = brain_extraction( input_image ) 1654 1655 tonlycerebrum = brain_extraction( input_template ) 1656 template = ants.rank_intensity( input_template ) 1657 1658 regsegits=[200,200,20] 1659 1660# # upsample the template if we are passing SR as input 1661# if min(ants.get_spacing(img)) < 0.8 : 1662# regsegits=[200,200,200,20] 1663# template = ants.resample_image( template, (0.5,0.5,0.5), interp_type = 0 ) 1664# tonlycerebrum = ants.resample_image_to_target( tonlycerebrum, 1665# template, 1666# interp_type='nearestNeighbor', 1667# ) 1668 1669 if is_test: 1670 regsegits=[8,0,0] 1671 1672 input_template_hemisphere_labels = ants.resample_image_to_target( 1673 input_template_hemisphere_labels, 1674 template, 1675 interp_type='nearestNeighbor', 1676 ) 1677 1678 # now do a hemisphere focused registration 1679 synL = localsyn( 1680 img=img*ionlycerebrum, 1681 template=template*tonlycerebrum, 1682 hemiS=input_image_hemisphere_segmentation, 1683 templateHemi=input_template_hemisphere_labels, 1684 whichHemi=1, 1685 padder=padding, 1686 iterations=regsegits, 1687 output_prefix = output_prefix + "left_hemi_reg", 1688 total_sigma=total_sigma, 1689 ) 1690 synR = localsyn( 1691 img=img*ionlycerebrum, 1692 template=template*tonlycerebrum, 1693 hemiS=input_image_hemisphere_segmentation, 1694 templateHemi=input_template_hemisphere_labels, 1695 whichHemi=2, 1696 padder=padding, 1697 iterations=regsegits, 1698 output_prefix = output_prefix + "right_hemi_reg", 1699 total_sigma=total_sigma, 1700 ) 1701 1702 ants.image_write(synL['warpedmovout'], output_prefix + "left_hemi_reg.nii.gz" ) 1703 ants.image_write(synR['warpedmovout'], output_prefix + "right_hemi_reg.nii.gz" ) 1704 1705 fignameL = output_prefix + "_left_hemi_reg.png" 1706 ants.plot(synL['warpedmovout'],axis=2,ncol=8,nslices=24,filename=fignameL, black_bg=False, crop=True ) 1707 1708 fignameR = output_prefix + "_right_hemi_reg.png" 1709 ants.plot(synR['warpedmovout'],axis=2,ncol=8,nslices=24,filename=fignameR, black_bg=False, crop=True ) 1710 1711 lhjac = ants.create_jacobian_determinant_image( 1712 synL['warpedmovout'], 1713 synL['fwdtransforms'][0], 1714 do_log=1 1715 ) 1716 ants.image_write( lhjac, output_prefix+'left_hemi_jacobian.nii.gz' ) 1717 1718 rhjac = ants.create_jacobian_determinant_image( 1719 synR['warpedmovout'], 1720 synR['fwdtransforms'][0], 1721 do_log=1 1722 ) 1723 ants.image_write( rhjac, output_prefix+'right_hemi_jacobian.nii.gz' ) 1724 return { 1725 "synL":synL, 1726 "synLpng":fignameL, 1727 "synR":synR, 1728 "synRpng":fignameR, 1729 "lhjac":lhjac, 1730 "rhjac":rhjac 1731 }
hemisphere focused registration that will produce jacobians and figures to support data inspection
input_image: input image
input_image_tissue_segmentation: segmentation produced in ANTs style ie with labels defined by atropos brain segmentation (1 to 6)
input_image_hemisphere_segmentation: left (1) and right (2) hemisphere segmentation
input_template: template to which we register; prefer a population-specific relatively high-resolution template instead of MNI or biobank.
input_template_hemisphere_labels: a segmentation image of template hemispheres with left labeled 1 and right labeled 2
output_prefix: a path and prefix for registration related outputs
padding: number of voxels to pad images, needed for diffzero
labels_to_register: list of integer segmentation labels to use to define the tissue types / regions of the brain to register.
total_sigma: scalar >= 0.0 ; higher means more constrained registration.
is_test: boolean. this function can be long running by default. this would help testing more quickly by running fewer iterations.
def
region_reg( input_image, input_image_tissue_segmentation, input_image_region_segmentation, input_template, input_template_region_segmentation, output_prefix, padding=10, total_sigma=0.5, is_test=False):
1735def region_reg( 1736 input_image, 1737 input_image_tissue_segmentation, 1738 input_image_region_segmentation, 1739 input_template, 1740 input_template_region_segmentation, 1741 output_prefix, 1742 padding=10, 1743 total_sigma=0.5, 1744 is_test=False ): 1745 """ 1746 region focused registration that will produce jacobians and figures to 1747 support data inspection. region-defining images should be binary. 1748 1749 input_image: input image 1750 1751 input_image_tissue_segmentation: segmentation produced in ANTs style ie with 1752 labels defined by atropos brain segmentation (1 to 6) 1753 1754 input_image_region_segmentation: a local region to register - binary. 1755 1756 input_template: template to which we register; prefer a population-specific 1757 relatively high-resolution template instead of MNI or biobank. brain extracted. 1758 1759 input_template_region_segmentation: a segmentation image of template regions - binary. 1760 1761 output_prefix: a path and prefix for registration related outputs 1762 1763 padding: number of voxels to pad images, needed for diffzero 1764 1765 total_sigma: scalar >= 0.0 ; higher means more constrained registration. 1766 1767 is_test: boolean. this function can be long running by default. this would 1768 help testing more quickly by running fewer iterations. 1769 1770 """ 1771 1772 img = ants.rank_intensity( input_image ) 1773 ionlycerebrum = brain_extraction( input_image ) 1774 template = ants.rank_intensity( input_template ) 1775 regsegits=[200,200,20] 1776 1777 # upsample the template if we are passing SR as input 1778 if min(ants.get_spacing(img)) < 0.8: 1779 regsegits=[200,200,200,20] 1780 template = ants.resample_image( template, (0.5,0.5,0.5), interp_type = 0 ) 1781 1782 if is_test: 1783 regsegits=[20,5,0] 1784 1785 input_template_region_segmentation = ants.resample_image_to_target( 1786 input_template_region_segmentation, 1787 template, 1788 interp_type='nearestNeighbor', 1789 ) 1790 1791 # now do a region focused registration 1792 synL = localsyn( 1793 img=img*ionlycerebrum, 1794 template=template, 1795 hemiS=input_image_region_segmentation, 1796 templateHemi=input_template_region_segmentation, 1797 whichHemi=1, 1798 padder=padding, 1799 iterations=regsegits, 1800 output_prefix = output_prefix + "region_reg", 1801 total_sigma=total_sigma, 1802 ) 1803 1804 ants.image_write(synL['warpedmovout'], output_prefix + "region_reg.nii.gz" ) 1805 1806 fignameL = output_prefix + "_region_reg.png" 1807 ants.plot(synL['warpedmovout'],axis=2,ncol=8,nslices=24,filename=fignameL, black_bg=False, crop=True ) 1808 1809 lhjac = ants.create_jacobian_determinant_image( 1810 synL['warpedmovout'], 1811 synL['fwdtransforms'][0], 1812 do_log=1 1813 ) 1814 ants.image_write( lhjac, output_prefix+'region_jacobian.nii.gz' ) 1815 1816 return { 1817 "synL":synL, 1818 "synLpng":fignameL, 1819 "lhjac":lhjac, 1820 'rankimg':img*ionlycerebrum, 1821 'template':template 1822 }
region focused registration that will produce jacobians and figures to support data inspection. region-defining images should be binary.
input_image: input image
input_image_tissue_segmentation: segmentation produced in ANTs style ie with labels defined by atropos brain segmentation (1 to 6)
input_image_region_segmentation: a local region to register - binary.
input_template: template to which we register; prefer a population-specific relatively high-resolution template instead of MNI or biobank. brain extracted.
input_template_region_segmentation: a segmentation image of template regions - binary.
output_prefix: a path and prefix for registration related outputs
padding: number of voxels to pad images, needed for diffzero
total_sigma: scalar >= 0.0 ; higher means more constrained registration.
is_test: boolean. this function can be long running by default. this would help testing more quickly by running fewer iterations.
def
t1_hypointensity( x, xsegmentation, xWMProbability, template, templateWMPrior, wmh_thresh=0.1):
1825def t1_hypointensity( x, xsegmentation, xWMProbability, template, templateWMPrior, wmh_thresh=0.1 ): 1826 """ 1827 provide measurements that may help decide if a given t1 image is likely 1828 to have hypointensity. 1829 1830 x: input image; bias-corrected, brain-extracted and denoised 1831 1832 xsegmentation: input image hard-segmentation results 1833 1834 xWMProbability: input image WM probability 1835 1836 template: template image 1837 1838 templateWMPrior: template tissue prior 1839 1840 wmh_thresh: float used to threshold WMH probability and produce summary data 1841 1842 returns: 1843 1844 - wmh_summary: summary data frame based on thresholding WMH probability at wmh_thresh 1845 - wmh_probability_image: probability image denoting WMH probability; higher values indicate 1846 that WMH is more likely 1847 - wmh_evidence_of_existence: an integral evidence that indicates the likelihood that the input 1848 image content supports the presence of white matter hypointensity. 1849 greater than zero is supportive of WMH. the higher, the more so. 1850 less than zero is evidence against. 1851 - wmh_max_prob: max probability of wmh 1852 - features: the features driving WMH predictons 1853 1854 """ 1855 1856 if False: # Need to retrain this model with refined inference code 1857 return { 1858 "wmh_summary":None, 1859 "wmh_probability_image":None, 1860 "wmh_evidence_of_existence":None, 1861 "wmh_max_prob":None, 1862 "features":None } 1863 1864 mybig = [88,128,128] 1865 templatesmall = ants.resample_image( template, mybig, use_voxels=True ) 1866 qaff = ants.registration( 1867 ants.rank_intensity(x), 1868 ants.rank_intensity(templatesmall), 'antsRegistrationSyNQuickRepro[s]', 1869 syn_sampling=2, 1870 syn_metric='CC', 1871 reg_iterations = [25,15,0,0],total_sigma=0.5, 1872 aff_metric='GC', random_seed=1 ) 1873 afftx = qaff['fwdtransforms'][1] 1874 templateWMPrior2x = ants.apply_transforms( x, templateWMPrior, qaff['fwdtransforms'] ) 1875 cerebrum = ants.threshold_image( xsegmentation, 2, 4 ) 1876 realWM = ants.threshold_image( templateWMPrior2x , 0.1, math.inf ) 1877 inimg = ants.rank_intensity( x ) 1878 parcellateWMdnz = ants.kmeans_segmentation( inimg, 2, realWM, mrf=0.3 )['probabilityimages'][0] 1879 x2template = ants.apply_transforms( templatesmall, inimg, afftx, whichtoinvert=[True] ) 1880 parcellateWMdnz2template = ants.apply_transforms( templatesmall, 1881 cerebrum * parcellateWMdnz, afftx, whichtoinvert=[True] ) 1882 # features = rank+dnz-image, lprob, wprob, wprior at mybig resolution 1883 f1 = x2template.numpy() 1884 f2 = parcellateWMdnz2template.numpy() 1885 f3 = ants.apply_transforms( templatesmall, xWMProbability, afftx, whichtoinvert=[True] ).numpy() 1886 f4 = ants.apply_transforms( templatesmall, templateWMPrior, qaff['fwdtransforms'][0] ).numpy() 1887 myfeatures = np.stack( (f1,f2,f3,f4), axis=3 ) 1888 newshape = np.concatenate( [ [1],np.asarray( myfeatures.shape )] ) 1889 myfeatures = myfeatures.reshape( newshape ) 1890 1891 inshape = [None,None,None,4] 1892 wmhunet = antspynet.create_unet_model_3d( inshape, 1893 number_of_outputs = 1, 1894 number_of_layers = 4, 1895 mode = 'sigmoid' ) 1896 1897 wmhunet.load_weights( get_data("simwmhseg", target_extension='.h5') ) 1898 1899 pp = wmhunet.predict( myfeatures ) 1900 1901 limg = ants.from_numpy( tf.squeeze( pp[0] ).numpy( ) ) 1902 limg = ants.copy_image_info( templatesmall, limg ) 1903 lesresam = ants.apply_transforms( x, limg, afftx, whichtoinvert=[False] ) 1904 # lesresam = lesresam * cerebrum 1905 rnmdl = antspynet.create_resnet_model_3d( inshape, 1906 number_of_outputs = 1, 1907 layers = (1,2,3), 1908 residual_block_schedule = (3,4,6,3), squeeze_and_excite = True, 1909 lowest_resolution = 32, cardinality = 1, mode = "regression" ) 1910 rnmdl.load_weights( get_data("simwmdisc", target_extension='.h5' ) ) 1911 qq = rnmdl.predict( myfeatures ) 1912 1913 lesresamb = ants.threshold_image( lesresam, wmh_thresh, 1.0 ) 1914 lgo=ants.label_geometry_measures( lesresamb, lesresam ) 1915 wmhsummary = pd.read_csv( get_data("wmh_evidence", target_extension='.csv' ) ) 1916 wmhsummary.at[0,'Value']=lgo.at[0,'VolumeInMillimeters'] 1917 wmhsummary.at[1,'Value']=lgo.at[0,'IntegratedIntensity'] 1918 wmhsummary.at[2,'Value']=float(qq) 1919 1920 return { 1921 "wmh_summary":wmhsummary, 1922 "wmh_probability_image":lesresam, 1923 "wmh_evidence_of_existence":float(qq), 1924 "wmh_max_prob":lesresam.max(), 1925 "features":myfeatures }
provide measurements that may help decide if a given t1 image is likely to have hypointensity.
x: input image; bias-corrected, brain-extracted and denoised
xsegmentation: input image hard-segmentation results
xWMProbability: input image WM probability
template: template image
templateWMPrior: template tissue prior
wmh_thresh: float used to threshold WMH probability and produce summary data
returns:
- wmh_summary: summary data frame based on thresholding WMH probability at wmh_thresh
- wmh_probability_image: probability image denoting WMH probability; higher values indicate
that WMH is more likely
- wmh_evidence_of_existence: an integral evidence that indicates the likelihood that the input
image content supports the presence of white matter hypointensity.
greater than zero is supportive of WMH. the higher, the more so.
less than zero is evidence against.
- wmh_max_prob: max probability of wmh
- features: the features driving WMH predictons
def
zoom_syn( target_image, template, template_segmentations, initial_registration, dilation=4, regIterations=[25]):
2996def zoom_syn( target_image, template, template_segmentations, 2997 initial_registration, 2998 dilation = 4, 2999 regIterations = [25] ): 3000 """ 3001 zoomed in syn - a hierarchical registration applied to a hierarchical segmentation 3002 3003 Initial registration is followed up by a refined and focused high-resolution registration. 3004 This is performed on the cropped image where the cropping region is determined 3005 by the first segmentation in the template_segmentations list. Segmentations 3006 after the first one are assumed to exist as sub-regions of the first. All 3007 segmentations are assumed to be binary. 3008 3009 Arguments 3010 --------- 3011 target_image : ants image at original resolution 3012 3013 template : ants image template to be mapped to the target image 3014 3015 template_segmentations : list of binary segmentation images 3016 3017 dilation : morphological dilation amount applied to the first segmentation and used for cropping 3018 3019 regIterations : parameter passed to ants.registration 3020 3021 Returns 3022 ------- 3023 dictionary 3024 containing segmentation and registration results in addition to cropping results 3025 3026 Example 3027 ------- 3028 >>> import ants 3029 >>> ireg = ants.registration( target_image, template, "antsRegistrationSyNQuickRepro[s]" ) 3030 >>> xxx = antspyt1w.zoom_syn( orb, template, level2segs, ireg ) 3031 """ 3032 croppertem = ants.iMath( template_segmentations[0], "MD", dilation ) 3033 templatecrop = ants.crop_image( template, croppertem ) 3034 cropper = ants.apply_transforms( target_image, 3035 croppertem, initial_registration['fwdtransforms'], 3036 interpolator='linear' ).threshold_image(0.5,1.e9) 3037 croplow = ants.crop_image( target_image, cropper ) 3038 synnerlow = ants.registration( croplow, templatecrop, 3039 'SyNOnly', gradStep = 0.20, regIterations = regIterations, randomSeed=1, 3040 syn_metric='cc', syn_sampling=2, total_sigma=0.5, 3041 initialTransform = initial_registration['fwdtransforms'] ) 3042 orlist = [] 3043 for jj in range(len(template_segmentations)): 3044 target_imageg = ants.apply_transforms( target_image, template_segmentations[jj], 3045 synnerlow['fwdtransforms'], 3046 interpolator='linear' ).threshold_image(0.5,1e9) 3047 orlist.append( target_imageg ) 3048 return{ 3049 'segmentations': orlist, 3050 'registration': synnerlow, 3051 'croppedimage': croplow, 3052 'croppingmask': cropper 3053 }
zoomed in syn - a hierarchical registration applied to a hierarchical segmentation
Initial registration is followed up by a refined and focused high-resolution registration. This is performed on the cropped image where the cropping region is determined by the first segmentation in the template_segmentations list. Segmentations after the first one are assumed to exist as sub-regions of the first. All segmentations are assumed to be binary.
Arguments
target_image : ants image at original resolution
template : ants image template to be mapped to the target image
template_segmentations : list of binary segmentation images
dilation : morphological dilation amount applied to the first segmentation and used for cropping
regIterations : parameter passed to ants.registration
Returns
dictionary containing segmentation and registration results in addition to cropping results
Example
>>> import ants
>>> ireg = ants.registration( target_image, template, "antsRegistrationSyNQuickRepro[s]" )
>>> xxx = antspyt1w.zoom_syn( orb, template, level2segs, ireg )
def
merge_hierarchical_csvs_to_wide_format( hierarchical_dataframes, col_names=None, identifier=None, identifier_name='u_hier_id', verbose=False):
3204def merge_hierarchical_csvs_to_wide_format( hierarchical_dataframes, col_names = None , identifier=None, identifier_name='u_hier_id', verbose=False ): 3205 """ 3206 standardized merging of output for dataframes produced by hierarchical function. 3207 3208 Arguments 3209 --------- 3210 hierarchical_dataframes : output of antspyt1w.hierarchical 3211 3212 identifier : unique subject identifier e.g. subject_001 3213 3214 identifier_name : string name for the unique identifer column e.g. subject_id 3215 3216 Returns 3217 ------- 3218 data frame in wide format 3219 3220 """ 3221 if identifier is None: 3222 identifier='A' 3223 wide_df = pd.DataFrame( ) 3224 icvkey='icv' 3225 if icvkey in hierarchical_dataframes.keys(): 3226 temp = hierarchical_dataframes[icvkey].copy() 3227 temp = temp.loc[:, ~temp.columns.str.contains('^Unnamed')] 3228 temp.insert(loc=0, column=identifier_name, value=identifier) 3229 temp = temp.set_index(identifier_name) 3230 wide_df = wide_df.join(temp,how='outer') 3231 for myvar in hierarchical_dataframes.keys(): 3232 if verbose: 3233 print( myvar ) 3234 if hierarchical_dataframes[myvar] is not None: 3235 jdf = hierarchical_dataframes[myvar].dropna(axis=0) 3236 jdf = jdf.loc[:, ~jdf.columns.str.contains('^Unnamed')] 3237 if col_names is not None : 3238 for col_name in col_names : 3239 if jdf.shape[0] > 1 and any( jdf.columns.str.contains(col_name)): 3240 varsofinterest = ["Description", col_name] 3241 jdfsub = jdf[varsofinterest] 3242 jdfsub.insert(loc=0, column=identifier_name, value=identifier) 3243 jdfsub = jdfsub.set_index([identifier_name, 'Description'])[col_name].unstack().add_prefix(col_name + '_') 3244 jdfsub.columns=jdfsub.columns 3245 jdfsub = jdfsub.rename(mapper=lambda x: x.strip().replace(' ', '_').lower(), axis=1) 3246 wide_df = wide_df.join(jdfsub,how='outer') 3247 if jdf.shape[0] > 1 and any( jdf.columns.str.contains('VolumeInMillimeters')): 3248 varsofinterest = ["Description", "VolumeInMillimeters"] 3249 jdfsub = jdf[varsofinterest] 3250 jdfsub.insert(loc=0, column=identifier_name, value=identifier) 3251 jdfsub=jdfsub.set_index([identifier_name, 'Description']).VolumeInMillimeters.unstack().add_prefix('Vol_') 3252 jdfsub.columns=jdfsub.columns+myvar 3253 jdfsub = jdfsub.rename(mapper=lambda x: x.strip().replace(' ', '_').lower(), axis=1) 3254 wide_df = wide_df.join(jdfsub,how='outer') 3255 if jdf.shape[0] > 1 and any( jdf.columns.str.contains('SurfaceAreaInMillimetersSquared')): 3256 varsofinterest = ["Description", "SurfaceAreaInMillimetersSquared"] 3257 jdfsub = jdf[varsofinterest] 3258 jdfsub.insert(loc=0, column=identifier_name, value=identifier) 3259 jdfsub=jdfsub.set_index([identifier_name, 'Description']).SurfaceAreaInMillimetersSquared.unstack().add_prefix('Area_') 3260 jdfsub.columns=jdfsub.columns+myvar 3261 jdfsub = jdfsub.rename(mapper=lambda x: x.strip().replace(' ', '_').lower(), axis=1) 3262 wide_df = wide_df.join(jdfsub,how='outer') 3263 if jdf.shape[0] > 1 and any( jdf.columns.str.contains('SurfaceAreaInMillimetersSquared')) and any( jdf.columns.str.contains('VolumeInMillimeters')): 3264 varsofinterest = ["Description", "VolumeInMillimeters", "SurfaceAreaInMillimetersSquared"] 3265 jdfsub = jdf[varsofinterest] 3266 jdfsub.insert(loc=0, column=identifier_name, value=identifier) 3267 jdfsub.insert(loc=1, column='thickness',value=jdfsub['VolumeInMillimeters']/jdfsub['SurfaceAreaInMillimetersSquared']) 3268 jdfsub=jdfsub.set_index([identifier_name, 'Description']).thickness.unstack().add_prefix('Thk_') 3269 jdfsub.columns=jdfsub.columns+myvar 3270 jdfsub = jdfsub.rename(mapper=lambda x: x.strip().replace(' ', '_').lower(), axis=1) 3271 wide_df = wide_df.join(jdfsub,how='outer') 3272 rbpkey='rbp' 3273 if rbpkey in hierarchical_dataframes.keys(): 3274 temp = hierarchical_dataframes[rbpkey].copy() 3275 temp = temp.loc[:, ~temp.columns.str.contains('^Unnamed')] 3276 temp.insert(loc=0, column=identifier_name, value=identifier) 3277 temp = temp.set_index(identifier_name) 3278 wide_df = wide_df.join(temp,how='outer') 3279 # handle wmh 3280 wmhkey='wmh' 3281 if wmhkey in hierarchical_dataframes.keys(): 3282 df=hierarchical_dataframes[wmhkey].copy() 3283 df.insert(loc=0, column=identifier_name, value=identifier) 3284 df = df.set_index(identifier_name) 3285 wmhval = df.loc[ df.Description == 'Volume_of_WMH','Value'] 3286 wide_df.insert(loc = 0, column = 'wmh_vol', value =wmhval ) 3287 wmhval = df.loc[ df.Description == 'Integral_WMH_probability','Value'] 3288 wide_df.insert(loc = 0, column = 'wmh_integral_prob', value =wmhval ) 3289 wmhval = df.loc[ df.Description == 'Log_Evidence','Value'] 3290 wide_df.insert(loc = 0, column = 'wmh_log_evidence', value =wmhval ) 3291 wide_df['wmh_log_evidence']=wmhval 3292 wide_df.insert(loc = 0, column = identifier_name, value = identifier) 3293 return wide_df
standardized merging of output for dataframes produced by hierarchical function.
Arguments
hierarchical_dataframes : output of antspyt1w.hierarchical
identifier : unique subject identifier e.g. subject_001
identifier_name : string name for the unique identifer column e.g. subject_id
Returns
data frame in wide format
def
map_intensity_to_dataframe(segmentation_type, intensity_image, segmentation_image):
283def map_intensity_to_dataframe(segmentation_type, intensity_image, segmentation_image): 284 """ 285 Match intensity values within segmentation labels to an appropriate data frame. 286 287 Arguments 288 --------- 289 segmentation_type : string or pandas.DataFrame 290 If string: 291 Name of the segmentation type data frame to retrieve using get_data. 292 Options include data available via get_data (e.g., "lobes"). 293 If pandas.DataFrame: 294 A pre-existing DataFrame containing label information. 295 296 intensity_image : ants.ANTsImage 297 ANTsImage with intensity values to summarize. 298 299 segmentation_image : ants.ANTsImage 300 ANTsImage with labels (or mostly the same values) as expected by segmentation_type. 301 302 Returns 303 ------- 304 pandas.DataFrame 305 A DataFrame summarizing the intensity values within segmentation labels. 306 """ 307 if isinstance(segmentation_type, str): 308 # Retrieve the data frame based on the segmentation type string 309 mydf_fn = get_data(segmentation_type) 310 mydf = pd.read_csv(mydf_fn) 311 elif isinstance(segmentation_type, pd.DataFrame): 312 # Use the provided DataFrame 313 mydf = segmentation_type 314 else: 315 raise ValueError("segmentation_type must be either a string or a pandas DataFrame") 316 317 # Get label statistics for the intensity and segmentation images 318 mylgo = ants.label_stats(intensity_image, segmentation_image) 319 320 # Rename the label column to match expected name for merging 321 mylgo = mylgo.rename(columns={'LabelValue': 'Label'}) 322 323 # Merge the segmentation DataFrame with the label statistics 324 result = pd.merge(mydf, mylgo, how='left', on=["Label"]) 325 326 return result
Match intensity values within segmentation labels to an appropriate data frame.
Arguments
segmentation_type : string or pandas.DataFrame If string: Name of the segmentation type data frame to retrieve using get_data. Options include data available via get_data (e.g., "lobes"). If pandas.DataFrame: A pre-existing DataFrame containing label information.
intensity_image : ants.ANTsImage ANTsImage with intensity values to summarize.
segmentation_image : ants.ANTsImage ANTsImage with labels (or mostly the same values) as expected by segmentation_type.
Returns
pandas.DataFrame A DataFrame summarizing the intensity values within segmentation labels.
def
deep_nbm( t1, nbm_weights, binary_mask=None, deform=False, aged_template=False, csfquantile=None, reflect=False, verbose=False):
1929def deep_nbm( t1, 1930 nbm_weights, 1931 binary_mask=None, 1932 deform=False, 1933 aged_template=False, 1934 csfquantile=None, 1935 reflect=False, 1936 verbose=False ): 1937 1938 """ 1939 CH13 and Nucleus basalis of Meynert segmentation and subdivision 1940 1941 Perform CH13 and NBM segmentation in T1 images using Avants labels. 1942 1943 t1 : T1-weighted neuroimage antsImage - already brain extracted 1944 1945 nbm_weights : string weight file for parcellating unet 1946 1947 binary_mask : will restrict output to this mask 1948 1949 deform : boolean to correct for image deformation 1950 1951 aged_template : boolean to control which template to use 1952 1953 csfquantile: if not None, will try to remove CSF from the image. 1954 0.25 may be a good value to try. 1955 1956 verbose: boolean 1957 1958 The labeling is as follows: 1959 1960 Label,Description,Side 1961 1,CH13_left,left 1962 2,CH13_right,right 1963 3,NBM_left_ant,left 1964 4,NBM_left_mid,left 1965 5,NBM_left_pos,left 1966 6,NBM_right_ant,right 1967 7,NBM_right_mid,right 1968 8,NBM_right_pos,right 1969 1970 Failure modes will include odd image orientation (in which case you might 1971 use the registration option). A more nefarious issue can be a poor extraction 1972 of the cerebrum in the inferior frontal lobe. These can be unpredictable 1973 but if one sees a bad extraction, please check the mask that is output by 1974 this function to see if it excludes non-cerebral tissue. 1975 1976 """ 1977 1978 if not aged_template: 1979 refimg = ants.image_read( get_data( "CIT168_T1w_700um_pad", target_extension='.nii.gz' )) 1980 refimg = ants.rank_intensity( refimg ) 1981 refimg = ants.resample_image( refimg, [0.5,0.5,0.5] ) 1982 refimgseg = ants.image_read( get_data( "CIT168_basal_forebrain", target_extension='.nii.gz' )) 1983 refimgsmall = ants.resample_image( refimg, [2.0,2.0,2.0] ) 1984 else: 1985 refimg = ants.image_read( get_data( "CIT168_T1w_700um_pad_adni", target_extension='.nii.gz' )) 1986 refimg = ants.rank_intensity( refimg ) 1987 refimg = ants.resample_image( refimg, [0.5,0.5,0.5] ) 1988 refimgseg = ants.image_read( get_data( "CIT168_basal_forebrain_adni", target_extension='.nii.gz' )) 1989 refimgsmall = ants.resample_image( refimg, [2.0,2.0,2.0] ) 1990 1991 pt_labels = [1,2,3,4] 1992 group_labels = [0,1,2,3,4,5,6,7,8] 1993 reflection_labels = [0,2,1,6,7,8,3,4,5] 1994 crop_size = [144,96,64] 1995 1996 def nbmpreprocess( img, pt_labels, group_labels, masker=None, csfquantile=None, returndef=False, reflect=False ): 1997 1998 if masker is None: 1999 imgr = ants.rank_intensity( img ) 2000 else: 2001 imgr = ants.rank_intensity( img, mask = masker ) 2002 2003 if csfquantile is not None and masker is None: 2004 masker = ants.threshold_image( imgr, np.quantile(imgr[imgr>1e-4], csfquantile ), 1e9 ) 2005 2006 if masker is not None: 2007 imgr = imgr * masker 2008 2009 imgrsmall = ants.resample_image( imgr, [1,1,1] ) 2010 if not reflect: 2011 reg = ants.registration( 2012 refimgsmall, 2013 imgrsmall, 2014 'antsRegistrationSyNQuickRepro[s]', 2015 reg_iterations = [200,200,20],total_sigma=0.5, 2016 verbose=False ) 2017 else: 2018 myref = ants.reflect_image( imgrsmall, axis=0, tx='Translation' ) 2019 reg = ants.registration( 2020 refimgsmall, 2021 imgrsmall, 2022 'antsRegistrationSyNQuickRepro[s]', 2023 reg_iterations = [200,200,20], 2024 total_sigma=0.5, 2025 initial_transform = myref['fwdtransforms'][0], 2026 verbose=False ) 2027 if not returndef: 2028 imgraff = ants.apply_transforms( refimg, imgr, reg['fwdtransforms'][1], interpolator='linear' ) 2029 imgseg = ants.apply_transforms( refimg, refimgseg, reg['invtransforms'][1], interpolator='nearestNeighbor' ) 2030 else: 2031 imgraff = ants.apply_transforms( refimg, imgr, reg['fwdtransforms'], interpolator='linear' ) 2032 imgseg = ants.image_clone( refimgseg ) 2033 binseg = ants.mask_image( imgseg, imgseg, pt_labels, binarize=True ) 2034 imgseg = ants.mask_image( imgseg, imgseg, group_labels, binarize=False ) 2035 com = ants.get_center_of_mass( binseg ) 2036 return { 2037 "img": imgraff, 2038 "seg": imgseg, 2039 "imgc": special_crop( imgraff, com, crop_size ), 2040 "segc": special_crop( imgseg, com, crop_size ), 2041 "reg" : reg, 2042 "mask": masker 2043 } 2044 2045 2046 nLabels = len( group_labels ) 2047 number_of_outputs = len(group_labels) 2048 number_of_channels = 1 2049 ################################################ 2050 unet0 = antspynet.create_unet_model_3d( 2051 [ None, None, None, number_of_channels ], 2052 number_of_outputs = 1, # number of landmarks must be known 2053 number_of_layers = 4, # should optimize this wrt criterion 2054 number_of_filters_at_base_layer = 32, # should optimize this wrt criterion 2055 convolution_kernel_size = 3, # maybe should optimize this wrt criterion 2056 deconvolution_kernel_size = 2, 2057 pool_size = 2, 2058 strides = 2, 2059 dropout_rate = 0.0, 2060 weight_decay = 0, 2061 additional_options = "nnUnetActivationStyle", 2062 mode = "sigmoid" ) 2063 2064 unet1 = antspynet.create_unet_model_3d( 2065 [None,None,None,2], 2066 number_of_outputs=number_of_outputs, 2067 mode="classification", 2068 number_of_filters=(32, 64, 96, 128, 256), 2069 convolution_kernel_size=(3, 3, 3), 2070 deconvolution_kernel_size=(2, 2, 2), 2071 dropout_rate=0.0, 2072 weight_decay=0, 2073 additional_options = "nnUnetActivationStyle") 2074 2075 # concat output to input and pass to 2nd net 2076 # nextin = tf.concat( [ unet0.inputs[0], unet0.outputs[0] ], axis=4 ) 2077 import tensorflow as tf 2078 from tensorflow.keras.layers import Layer 2079 class myConcat(Layer): 2080 def call(self, x): 2081 return tf.concat(x, axis=4 ) 2082 nextin = myConcat()( [ unet0.inputs[0], unet0.outputs[0] ] ) 2083 2084 unetonnet = unet1( nextin ) 2085 unet_model = tf.keras.models.Model( 2086 unet0.inputs, 2087 [ unetonnet, unet0.outputs[0] ] ) 2088 2089 unet_model.load_weights( nbm_weights ) 2090 2091 imgprepro = nbmpreprocess( t1, pt_labels, group_labels, 2092 csfquantile = csfquantile, 2093 returndef = deform, 2094 reflect = reflect ) 2095 2096 def map_back( relo, t1, imgprepro, interpolator='linear', deform = False ): 2097 if not deform: 2098 relo = ants.apply_transforms( t1, relo, 2099 imgprepro['reg']['invtransforms'][0], whichtoinvert=[True], 2100 interpolator=interpolator ) 2101 else: 2102 relo = ants.apply_transforms( t1, relo, 2103 imgprepro['reg']['invtransforms'], interpolator=interpolator ) 2104 return relo 2105 2106 2107 #################################################### 2108 physspaceBF = imgprepro['imgc'] 2109 if verbose: 2110 print( "Mean preprocessed cropped image in nbm seg " + str(physspaceBF.mean() ) ) 2111 2112 tfarr1 = tf.cast( physspaceBF.numpy() ,'float32' ) 2113 newshapeBF = list( tfarr1.shape ) 2114 newshapeBF.insert(0,1) 2115 newshapeBF.insert(4,1) 2116 tfarr1 = tf.reshape(tfarr1, newshapeBF ) 2117 snpred = unet_model.predict( tfarr1 ) 2118 segpred = snpred[0] 2119 sigmoidpred = snpred[1] 2120 snpred1_image = ants.from_numpy( sigmoidpred[0,:,:,:,0] ) 2121 snpred1_image = ants.copy_image_info( physspaceBF, snpred1_image ) 2122 snpred1_image = map_back( snpred1_image, t1, imgprepro, 'linear', deform ) 2123 bint = ants.threshold_image( snpred1_image, 0.5, 1.0 ) 2124 probability_images = [] 2125 for jj in range(number_of_outputs-1): 2126 temp = ants.from_numpy( segpred[0,:,:,:,jj+1] ) 2127 temp = ants.copy_image_info( physspaceBF, temp ) 2128 temp = map_back( temp, t1, imgprepro, 'linear', deform ) 2129 probability_images.append( temp ) 2130 image_matrix = ants.image_list_to_matrix(probability_images, bint) 2131 segmentation_matrix = (np.argmax(image_matrix, axis=0) + 1) 2132 segmentation_image = ants.matrix_to_images(np.expand_dims(segmentation_matrix, axis=0), bint)[0] 2133 relabeled_image = ants.image_clone(segmentation_image) 2134 if not reflect: 2135 for i in range(1,len(group_labels)): 2136 relabeled_image[segmentation_image==(i)] = group_labels[i] 2137 else: 2138 for i in range(1,len(group_labels)): 2139 relabeled_image[segmentation_image==(i)] = reflection_labels[i] 2140 2141 bfsegdesc = map_segmentation_to_dataframe( 'nbm3CH13', relabeled_image ) 2142 2143 return { 'segmentation':relabeled_image, 2144 'description':bfsegdesc, 2145 'mask': imgprepro['mask'], 2146 'probability_images': probability_images }
CH13 and Nucleus basalis of Meynert segmentation and subdivision
Perform CH13 and NBM segmentation in T1 images using Avants labels.
t1 : T1-weighted neuroimage antsImage - already brain extracted
nbm_weights : string weight file for parcellating unet
binary_mask : will restrict output to this mask
deform : boolean to correct for image deformation
aged_template : boolean to control which template to use
csfquantile: if not None, will try to remove CSF from the image. 0.25 may be a good value to try.
verbose: boolean
The labeling is as follows:
Label,Description,Side 1,CH13_left,left 2,CH13_right,right 3,NBM_left_ant,left 4,NBM_left_mid,left 5,NBM_left_pos,left 6,NBM_right_ant,right 7,NBM_right_mid,right 8,NBM_right_pos,right
Failure modes will include odd image orientation (in which case you might use the registration option). A more nefarious issue can be a poor extraction of the cerebrum in the inferior frontal lobe. These can be unpredictable but if one sees a bad extraction, please check the mask that is output by this function to see if it excludes non-cerebral tissue.
map_cit168
def
deep_cit168( t1, binary_mask=None, syn_type='antsRegistrationSyNQuickRepro[s]', priors=None, verbose=False):
2385def deep_cit168( t1, binary_mask = None, 2386 syn_type='antsRegistrationSyNQuickRepro[s]', 2387 priors = None, verbose = False): 2388 2389 """ 2390 CIT168 atlas segmentation with a parcellation unet. 2391 2392 Perform CIT168 segmentation in T1 images using Pauli atlas (CIT168) labels. 2393 2394 t1 : T1-weighted neuroimage antsImage - already brain extracted. image should 2395 be normalized 0 to 1 and with a nicely uniform histogram (no major outliers). 2396 we do a little work internally to help this but no guarantees it will handle 2397 all possible confounding issues. 2398 2399 binary_mask : will restrict output to this mask 2400 2401 syn_type : the type of registration used for generating priors; usually 2402 either SyN or antsRegistrationSyNQuickRepro[s] for repeatable results 2403 2404 priors : the user can provide their own priors through this argument; for 2405 example, the user may run this function twice, with the output of the first 2406 giving input to the second run. 2407 2408 verbose: boolean 2409 2410 Failure modes will primarily occur around red nucleus and caudate nucleus. 2411 For the latter, one might consider masking by the ventricular CSF, in particular 2412 near the anterior and inferior portion of the caudate in subjects with 2413 large ventricles. Low quality images with high atropy are also likely outside 2414 of the current range of the trained models. Iterating the model may help. 2415 2416 """ 2417 def tfsubset( x, indices ): 2418 with tf.device('/CPU:0'): 2419 outlist = [] 2420 for k in indices: 2421 outlist.append( x[:,:,:,int(k)] ) 2422 return tf.stack( outlist, axis=3 ) 2423 2424 def tfsubsetbatch( x, indices ): 2425 with tf.device('/CPU:0'): 2426 outlist2 = [] 2427 for j in range( len( x ) ): 2428 outlist = [] 2429 for k in indices: 2430 if len( x[j].shape ) == 5: 2431 outlist.append( x[j][k,:,:,:,:] ) 2432 if len( x[j].shape ) == 4: 2433 outlist.append( x[j][k,:,:,:] ) 2434 outlist2.append( tf.stack( outlist, axis=0 ) ) 2435 return outlist2 2436 2437 2438 registration = True 2439 cit168seg = t1 * 0 2440 myprior = ants.image_read(get_data("det_atlas_25_pad_LR_adni", target_extension=".nii.gz")) 2441 nbmtemplate = ants.image_read( get_data( "CIT168_T1w_700um_pad_adni", target_extension=".nii.gz" ) ) 2442 nbmtemplate = ants.resample_image( nbmtemplate, [0.5,0.5,0.5] ) 2443 templateSmall = ants.resample_image( nbmtemplate, [2.0,2.0,2.0] ) 2444 orireg = ants.registration( 2445 fixed = templateSmall, 2446 moving = ants.iMath( t1, "Normalize" ), 2447 type_of_transform=syn_type, total_sigma=0.5, verbose=False ) 2448 image = ants.apply_transforms( nbmtemplate, ants.iMath( t1, "Normalize" ), 2449 orireg['fwdtransforms'][1] ) 2450 image = ants.iMath( image, "TruncateIntensity",0.001,0.999).iMath("Normalize") 2451 patchSize = [ 160,160,112 ] 2452 if priors is None: 2453 priortosub = ants.apply_transforms( image, myprior, 2454 orireg['invtransforms'][1], interpolator='nearestNeighbor' ) 2455 else: 2456 if verbose: 2457 print("using priors") 2458 priortosub = ants.apply_transforms( image, priors, 2459 orireg['fwdtransforms'][1], interpolator='nearestNeighbor' ) 2460 bmask = ants.threshold_image( priortosub, 1, 999 ) 2461 # this identifies the cropping location - assumes a good registration 2462 pt = list( ants.get_center_of_mass( bmask ) ) 2463 pt[1] = pt[1] + 10.0 # same as we did in training 2464 2465 physspaceCIT = special_crop( image, pt, patchSize) # image 2466 2467 if binary_mask is not None: 2468 binary_mask_use = ants.apply_transforms( nbmtemplate, binary_mask, 2469 orireg['fwdtransforms'][1] ) 2470 binary_mask_use = special_crop( binary_mask_use, pt, patchSize) 2471 2472 for sn in [True,False]: 2473 if sn: 2474 group_labels = [0,7,8,9,23,24,25,33,34] 2475 newfn=get_data( "deepCIT168_sn", target_extension=".h5" ) 2476 else: 2477 group_labels = [0,1,2,5,6,17,18,21,22] 2478 newfn=get_data( "deepCIT168", target_extension=".h5" ) 2479 2480 number_of_outputs = len(group_labels) 2481 number_of_channels = len(group_labels) 2482 2483 unet0 = antspynet.create_unet_model_3d( 2484 [ None, None, None, number_of_channels ], 2485 number_of_outputs = 1, # number of landmarks must be known 2486 number_of_layers = 4, # should optimize this wrt criterion 2487 number_of_filters_at_base_layer = 32, # should optimize this wrt criterion 2488 convolution_kernel_size = 3, # maybe should optimize this wrt criterion 2489 deconvolution_kernel_size = 2, 2490 pool_size = 2, 2491 strides = 2, 2492 dropout_rate = 0.0, 2493 weight_decay = 0, 2494 additional_options = "nnUnetActivationStyle", 2495 mode = "sigmoid" ) 2496 2497 unet1 = antspynet.create_unet_model_3d( 2498 [None,None,None,2], 2499 number_of_outputs=number_of_outputs, 2500 mode="classification", 2501 number_of_filters=(32, 64, 96, 128, 256), 2502 convolution_kernel_size=(3, 3, 3), 2503 deconvolution_kernel_size=(2, 2, 2), 2504 dropout_rate=0.0, 2505 weight_decay=0, 2506 additional_options = "nnUnetActivationStyle") 2507 2508# temp = tf.split( unet0.inputs[0], 9, axis=4 ) 2509 import tensorflow as tf 2510 from tensorflow.keras.layers import Layer 2511 class mySplit(Layer): 2512 def call(self, x): 2513 return tf.split(x, 9, axis=4 ) 2514 temp = mySplit()( unet0.inputs[0] ) 2515 2516 temp[1] = unet0.outputs[0] 2517 2518 class myConcat(Layer): 2519 def call(self, x): 2520 return tf.concat(x, axis=4 ) 2521# newmult = tf.concat( temp[0:2], axis=4 ) 2522 newmult = myConcat()( temp[0:2] ) 2523 unetonnet = unet1( newmult ) 2524 unet_model = tf.keras.models.Model( 2525 unet0.inputs, 2526 [ unetonnet, unet0.outputs[0] ] ) 2527 unet_model.load_weights( newfn ) 2528 ################### 2529 nbmprior = special_crop( priortosub, pt, patchSize).numpy() # prior 2530 imgnp = tf.reshape( physspaceCIT.numpy(), [160, 160, 112,1]) 2531 nbmprior = tf.one_hot( nbmprior, 35 ) 2532 nbmprior = tfsubset( nbmprior, group_labels[1:len(group_labels)] ) 2533 imgnp = tf.reshape( tf.concat( [imgnp,nbmprior], axis=3), [1,160, 160, 112,9]) 2534 2535 nbmpred = unet_model.predict( imgnp ) 2536 segpred = nbmpred[0] 2537 sigmoidpred = nbmpred[1] 2538 2539 def map_back_cit( relo, t1, orireg, interpolator='linear' ): 2540 relo = ants.apply_transforms( t1, relo, 2541 orireg['invtransforms'][0], whichtoinvert=[True], 2542 interpolator=interpolator ) 2543 return relo 2544 2545 nbmpred1_image = ants.from_numpy( sigmoidpred[0,:,:,:,0] ) 2546 nbmpred1_image = ants.copy_image_info( physspaceCIT, nbmpred1_image ) 2547 nbmpred1_image = map_back_cit( nbmpred1_image, t1, orireg, 'linear' ) 2548 if binary_mask is not None: 2549 nbmpred1_image = nbmpred1_image * binary_mask 2550 bint = ants.threshold_image( nbmpred1_image, 0.5, 1.0 ) 2551 2552 probability_images = [] 2553 for jj in range(1,len(group_labels)): 2554 temp = ants.from_numpy( segpred[0,:,:,:,jj] ) 2555 temp = ants.copy_image_info( physspaceCIT, temp ) 2556 temp = map_back_cit( temp, t1, orireg, 'linear' ) 2557 probability_images.append( temp ) 2558 2559 image_matrix = ants.image_list_to_matrix(probability_images, bint) 2560 segmentation_matrix = (np.argmax(image_matrix, axis=0) + 1) 2561 segmentation_image = ants.matrix_to_images(np.expand_dims(segmentation_matrix, axis=0), bint)[0] 2562 relabeled_image = ants.image_clone(segmentation_image*0.) 2563 for i in np.unique(segmentation_image.numpy()): 2564 if i > 0 : 2565 temp = ants.threshold_image(segmentation_image,i,i) 2566 if group_labels[int(i)] < 33: 2567 temp = ants.iMath( temp, "GetLargestComponent",1) 2568 relabeled_image = relabeled_image + temp*group_labels[int(i)] 2569 cit168seg = cit168seg + relabeled_image 2570 2571 cit168segdesc = map_segmentation_to_dataframe( 'CIT168_Reinf_Learn_v1_label_descriptions_pad', cit168seg ).dropna(axis=0) 2572 2573 return { 'segmentation':cit168seg, 'description':cit168segdesc }
CIT168 atlas segmentation with a parcellation unet.
Perform CIT168 segmentation in T1 images using Pauli atlas (CIT168) labels.
t1 : T1-weighted neuroimage antsImage - already brain extracted. image should be normalized 0 to 1 and with a nicely uniform histogram (no major outliers). we do a little work internally to help this but no guarantees it will handle all possible confounding issues.
binary_mask : will restrict output to this mask
syn_type : the type of registration used for generating priors; usually either SyN or antsRegistrationSyNQuickRepro[s] for repeatable results
priors : the user can provide their own priors through this argument; for example, the user may run this function twice, with the output of the first giving input to the second run.
verbose: boolean
Failure modes will primarily occur around red nucleus and caudate nucleus. For the latter, one might consider masking by the ventricular CSF, in particular near the anterior and inferior portion of the caudate in subjects with large ventricles. Low quality images with high atropy are also likely outside of the current range of the trained models. Iterating the model may help.
def
minimal_sr_preprocessing(x, imgbxt=None):
3657def minimal_sr_preprocessing( x, imgbxt=None ): 3658 """ 3659 x : input t1 image (raw) 3660 3661 imgbxt : optional existing brain extraction 3662 3663 outputs: preprocessedT1, hemisphereLabels 3664 """ 3665 if x.dimension != 3: 3666 raise ValueError('hierarchical: input image should be 3-dimensional') 3667 3668 tfn = get_data('T_template0', target_extension='.nii.gz' ) 3669 tlrfn = get_data('T_template0_LR', target_extension='.nii.gz' ) 3670 3671 ##### read images and do simple bxt ops 3672 templatea = ants.image_read( tfn ) 3673 templatea = ( templatea * antspynet.brain_extraction( templatea, 't1' ) ).iMath( "Normalize" ) 3674 templatealr = ants.image_read( tlrfn ) 3675 if imgbxt is None: 3676 imgbxt = brain_extraction( ants.iMath( x, "Normalize" ) ) 3677 img = preprocess_intensity( ants.iMath( x, "Normalize" ), imgbxt ) 3678 img = ants.iMath( img, "Normalize" ) 3679 mylr = label_hemispheres( img, templatea, templatealr ) 3680 return img, mylr * imgbxt
x : input t1 image (raw)
imgbxt : optional existing brain extraction
outputs: preprocessedT1, hemisphereLabels