In [1]: 
import SimpleITK as sitk
from myshow import myshow, myshow3d
In [2]: 
img_T1 = sitk.ReadImage("Data/nac-brain-atlas-1.0/volumes/A1_grayT1.nrrd")
img_T2 = sitk.ReadImage("Data/nac-brain-atlas-1.0/volumes/A1_grayT2.nrrd")
img_T1_255 = sitk.Cast(sitk.RescaleIntensity(img_T1), sitk.sitkUInt8)
img_T1_255 = sitk.Cast(sitk.RescaleIntensity(img_T1), sitk.sitkUInt8)
In [3]: 
sitk.Show(img_T1, title="T1")
In [4]: 
idx = (106,116,141)
pt = img_T1.TransformIndexToPhysicalPoint(idx)
In [5]: 
seg = sitk.Image(img_T1.GetSize(), sitk.sitkUInt8)
seg.CopyInformation(img_T1)
seg[idx] = 1
seg = sitk.BinaryDilate(seg, 3)
myshow3d(sitk.LabelOverlay(img_T1_255, seg), zslices=range(idx[2]-3, idx[2]+4, 3), dpi=30, title="Initial Seed")
In [6]: 
stats = sitk.LabelStatisticsImageFilter()
stats.Execute(img_T1, seg)
print stats

itk::simple::LabelStatisticsImageFilter
	 Label(0)[Count] = 1.6777e+07
	 Label(0)[Maximum] = 2126
	 Label(0)[Mean] = 87.8872
	 Label(0)[Minimum] = -122
	 Label(0)[Sigma] = 176.613
	 Label(0)[Sum] = 1.47449e+09
	 Label(0)[Variance] = 31192.1
	 Label(0)[approxMedian] = -3.45312
	 Label(0)[BoundingBox] = [ 0, 255, 0, 255, 0, 255 ]
	 Label(1)[Count] = 179
	 Label(1)[Maximum] = 444
	 Label(1)[Mean] = 382.369
	 Label(1)[Minimum] = 313
	 Label(1)[Sigma] = 18.1829
	 Label(1)[Sum] = 68444
	 Label(1)[Variance] = 330.616
	 Label(1)[approxMedian] = 382.922
	 Label(1)[BoundingBox] = [ 103, 109, 113, 119, 138, 144 ]
In [7]: 
factor = 1.5
lower_threshold = stats.GetMean(1)-factor*stats.GetSigma(1)
upper_threshold = stats.GetMean(1)+factor*stats.GetSigma(1)
In [8]: 
init_ls = sitk.SignedMaurerDistanceMap(seg, insideIsPositive=True, useImageSpacing=True)
In [9]: 
lsFilter = sitk.ThresholdSegmentationLevelSetImageFilter()
lsFilter.SetLowerThreshold(lower_threshold)
lsFilter.SetUpperThreshold(upper_threshold)
lsFilter.SetMaximumRMSError(0.02)
lsFilter.SetNumberOfIterations(100)
lsFilter.SetCurvatureScaling(1)
lsFilter.SetPropagationScaling(1)
lsFilter.ReverseExpansionDirectionOn()
Out[9]:
<SimpleITK.SimpleITK.ThresholdSegmentationLevelSetImageFilter; proxy of <Swig Object of type 'itk::simple::ThresholdSegmentationLevelSetImageFilter::Self *' at 0x11bab6210> >
In [10]: 
ls = lsFilter.Execute(init_ls, sitk.Cast(img_T1, sitk.sitkFloat32))
print lsFilter

itk::simple::ThresholdSegmentationLevelSetImageFilter
  LowerThreshold: 355.094
  UpperThreshold: 409.643
  MaximumRMSError: 0.02
  PropagationScaling: 1
  CurvatureScaling: 1
  NumberOfIterations: 100
  ReverseExpansionDirection: 1
  ElapsedIterations: 100
  RMSChange: 0.0369433
In [11]: 
zslice_offset = 4
t = "LevelSet after "+str(lsFilter.GetNumberOfIterations())+" iterations"
myshow3d(sitk.LabelOverlay(img_T1_255, ls > 0), zslices=range(idx[2]-zslice_offset,idx[2]+zslice_offset+1,zslice_offset), dpi=20, title=t)
In [12]: 
lsFilter.SetNumberOfIterations(25)
img_T1f = sitk.Cast(img_T1, sitk.sitkFloat32)
ls = init_ls
niter = 0
for i in range(0, 10):
    ls = lsFilter.Execute(ls, img_T1f)
    niter += lsFilter.GetNumberOfIterations()
    t = "LevelSet after "+str(niter)+" iterations and RMS "+str(lsFilter.GetRMSChange())
    fig = myshow3d(sitk.LabelOverlay(img_T1_255, ls > 0), zslices=range(idx[2]-zslice_offset,idx[2]+zslice_offset+1,zslice_offset), dpi=20, title=t)
In [12]: