Skip to content

Commit

Permalink
merge in Doug's changes from freesurfer/freesurfer/python/gems
Browse files Browse the repository at this point in the history 
  apply standalone samseg fixes
  • Loading branch information
@yhuang43
yhuang43 committed Jan 2, 2024
1 parent 3ca979f commit 60df589
Showing 1 changed file with 26 additions and 26 deletions.
52 changes: 26 additions & 26 deletions samseg/Samseg.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,16 +4,16 @@
import pickle
import scipy.io
import surfa as sf
from scipy.ndimage.morphology import binary_dilation as dilation
from scipy.ndimage import binary_dilation as dilation

import gems
from gems.utilities import Specification
from gems.BiasField import BiasField
from gems.ProbabilisticAtlas import ProbabilisticAtlas
from gems.GMM import GMM
from gems.Affine import Affine
from gems.SamsegUtility import *
from gems.merge_alphas import kvlMergeAlphas, kvlGetMergingFractionsTable
from samseg import gems
from .utilities import Specification
from .BiasField import BiasField
from .ProbabilisticAtlas import ProbabilisticAtlas
from .GMM import GMM
from .Affine import Affine
from .SamsegUtility import *
from .merge_alphas import kvlMergeAlphas, kvlGetMergingFractionsTable


eps = np.finfo(float).eps
Expand Down Expand Up @@ -74,10 +74,10 @@ def __init__(self,
raise ValueError('In photo mode, you cannot provide more than one input image volume')

input_vol = sf.load_volume(self.imageFileNames[0])
input_vol.save(self.savePath + '/original_input.mgz')
input_vol.save(os.path.join(self.savePath, 'original_input.mgz'))
if input_vol.nframes > 1:
input_vol = input_vol.mean(frames=True)
input_vol.save(self.savePath + '/grayscale_input.mgz')
input_vol.save(os.path.join(self.savePath, 'grayscale_input.mgz'))

# We also a small band of noise around the mask; otherwise the background/skull/etc may fit the cortex
self.photo_mask = input_vol.data > 0
Expand All @@ -87,26 +87,26 @@ def __init__(self,
rng = np.random.default_rng(2021)
input_vol.data[ring] = max_noise * rng.random(input_vol.data[ring].shape[0])
self.imageFileNames = []
self.imageFileNames.append(self.savePath + '/grayscale_input_with_ring.mgz')
self.imageFileNames.append(os.path.join(self.savePath,'grayscale_input_with_ring.mgz'))
input_vol.save(self.imageFileNames[0])


# Initialize some objects
self.affine = Affine( imageFileName=self.imageFileNames[0],
meshCollectionFileName=os.path.join(self.atlasDir, 'atlasForAffineRegistration.txt.gz'),
templateFileName=os.path.join(self.atlasDir, 'template.nii' ) )
templateFileName=os.path.join(self.atlasDir, 'template.nii.gz' ) )
self.probabilisticAtlas = ProbabilisticAtlas()

# Get full model specifications and optimization options (using default unless overridden by user)
# Note that, when processing photos, we point to a different GMM file by default!
self.optimizationOptions = getOptimizationOptions(atlasDir, userOptimizationOptions)
if dissectionPhoto and (gmmFileName is None):
if dissectionPhoto == 'left':
gmmFileName = self.atlasDir + '/photo.lh.sharedGMMParameters.txt'
gmmFileName = os.path.join(self.atlasDir, 'photo.lh.sharedGMMParameters.txt')
elif dissectionPhoto == 'right':
gmmFileName = self.atlasDir + '/photo.rh.sharedGMMParameters.txt'
gmmFileName = os.path.join(self.atlasDir, 'photo.rh.sharedGMMParameters.txt')
elif dissectionPhoto == 'both':
gmmFileName = self.atlasDir + '/photo.both.sharedGMMParameters.txt'
gmmFileName = os.path.join(self.atlasDir, 'photo.both.sharedGMMParameters.txt')
else:
sf.system.fatal('dissection photo mode must be left, right, or both')
self.modelSpecifications = getModelSpecifications(
Expand Down Expand Up @@ -229,16 +229,16 @@ def segment(self, costfile=None, timer=None, reg_only=False, transformFile=None,
if self.imageToImageTransformMatrix is None:

if self.dissectionPhoto is not None:
reference = self.savePath + '/grayscale_input.mgz'
reference = os.path.join(self.savePath, 'grayscale_input.mgz')
if self.dissectionPhoto=='left':
moving = self.atlasDir + '/exvivo.template.lh.suptent.nii'
moving = os.path.join(self.atlasDir, 'exvivo.template.lh.suptent.nii')
elif self.dissectionPhoto=='right':
moving = self.atlasDir + '/exvivo.template.rh.suptent.nii'
moving = os.path.join(self.atlasDir, 'exvivo.template.rh.suptent.nii')
elif self.dissectionPhoto=='both':
moving = self.atlasDir + '/exvivo.template.suptent.nii'
moving = os.path.join(self.atlasDir, 'exvivo.template.suptent.nii')
else:
sf.system.fatal('dissection photo mode must be left, right, or both')
transformFile = self.savePath + '/atlas2image.lta'
transformFile = os.path.join(self.savePath, 'atlas2image.lta')
cmd = 'mri_coreg --seed 2021 --mov ' + moving + ' --ref ' + reference + ' --reg ' + transformFile + \
' --dof 12 --threads ' + str(self.nthreads)
os.system(cmd)
Expand Down Expand Up @@ -318,7 +318,7 @@ def preProcess(self):
else:
self.imageBuffers, self.transform, self.voxelSpacing, self.cropping = readCroppedImages(
self.imageFileNames,
os.path.join(self.atlasDir, 'template.nii'),
os.path.join(self.atlasDir, 'template.nii.gz'),
self.imageToImageTransformMatrix
)

Expand Down Expand Up @@ -479,11 +479,11 @@ def writeResults(self, biasFields, posteriors):
print(self.scalingFactors[contrastNumber], file=fid)

else: # photos
self.writeImage(expBiasFields[..., 0], self.savePath + '/illlumination_field.mgz')
self.writeImage(expBiasFields[..., 0], os.path.join(self.savePath, 'illlumination_field.mgz'))
original_vol = sf.load_volume(self.originalImageFileNames[0])
bias_native = sf.load_volume(self.savePath + '/illlumination_field.mgz')
bias_native = sf.load_volume(os.path.join(self.savePath, 'illlumination_field.mgz'))
original_vol = original_vol / (1e-6 + bias_native)
original_vol.save(self.savePath + '/illlumination_corrected.mgz')
original_vol.save(os.path.join(self.savePath, 'illlumination_corrected.mgz'))

if self.savePosteriors:
posteriorPath = os.path.join(self.savePath, 'posteriors')
Expand Down Expand Up @@ -560,7 +560,7 @@ def saveWarpField(self, filename):

# extract geometries
source = sf.load_volume(self.imageFileNames[0]).geom
target = sf.load_volume(os.path.join(self.atlasDir, 'template.nii')).geom
target = sf.load_volume(os.path.join(self.atlasDir, 'template.nii.gz')).geom

# extract vox-to-vox template transform
# TODO: Grabbing the transform from the saved .mat file in either the cross or base
Expand Down

0 comments on commit 60df589

Please sign in to comment.