

In [1]:

    
import sys
sys.path.append('../code/functions')
sys.path.append('/home/simpleElastix/build/SimpleITK-build/Wrapping/Python')

import pickle
import cv2

import SimpleITK as sitk
import numpy as np
import matplotlib.pyplot as plt
import nibabel as nib

from tiffIO import loadTiff, unzipChannels
from connectLib import adaptiveThreshold

Load Data



In [2]:

    
tp1ChanList = unzipChannels(loadTiff('../data/SEP-GluA1-KI_tp1.tif'))



In [3]:

    
tp2ChanList = unzipChannels(loadTiff('../data/SEP-GluA1-KI_tp2.tif'))



In [ ]:

    
print tp1ChanList.shape
print tp2ChanList.shape



In [4]:

    
tp1Test = tp1ChanList[1][15:20]
tp2Test = tp2ChanList[1][15:20]



In [ ]:

    
print tp1Test.shape
print tp2Test.shape



In [ ]:

    
fig = plt.figure()
plt.imshow(tp1Test[0], cmap='gray')
plt.title('TP1 at z=10')
plt.show()



In [ ]:

    
fig = plt.figure()
plt.imshow(tp2Test[0], cmap='gray')
plt.title('TP2 at z=10')
plt.show()

Visualization Functions



In [5]:

    
def toDiff(imgA, imgB):
    ret = np.empty((imgA.shape[0], imgA.shape[1], 3), dtype=np.uint8)
    for y in range(imgA.shape[0]):
        for x in range(imgA.shape[1]):
            
            if imgA[y][x] and not imgB[y][x]:
                ret[y][x][0] = 255
                ret[y][x][1] = 0
                ret[y][x][2] = 0
            elif not imgA[y][x] and imgB[y][x]:
                ret[y][x][0] = 0
                ret[y][x][1] = 255
                ret[y][x][2] = 0
            elif imgA[y][x] and imgB[y][x]:
                ret[y][x][0] = 255
                ret[y][x][1] = 0
                ret[y][x][2] = 255
            else:
                ret[y][x][0] = 255
                ret[y][x][1] = 255
                ret[y][x][2] = 255
            
    return ret

def visDiff(sliceA, sliceB):
    disp = toDiff(sliceA, sliceB)
    return disp

def visVolDiff(volumeA, volumeB):
    for i in range(volumeA.shape[0]):
        plt.figure()
        plt.title('Disperity at z=' + str(i))
        plt.imshow(visDiff(volumeA[i], volumeB[i]))
        plt.show()

Strategy

First order of business is getting something that is an improvement on what we have currently. This should be simple; as what we have currently takes virtually forever to run.

I think the place to start is cleaning up the input images so the registration actually has something to hold on to



In [6]:

    
def preproc(img):
    binImg = adaptiveThreshold(img, 5, 5)
    return binImg



In [11]:

    
visVolDiff(tp1Test, tp2Test)



In [12]:

    
visVolDiff(preproc(tp1Test), preproc(tp2Test))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0

Now, that looks better. As a human I can now see how the images should pair, so (hopefully) this allows simple elastix to as well.

I plan to start simple - just a rigid transform with a low number of required valid samples to see where it gets us



In [37]:

    
%%time

SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['10']
pMap['RequiredRatioOfValidSamples'] = ['.05']

SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 26.7 s, sys: 412 ms, total: 27.1 s
Wall time: 5.75 s



In [32]:

    
transformParameterMap = SimpleElastix.GetTransformParameterMap()

transImg = sitk.GetImageFromArray(tp2Test)

imageFilter = sitk.ReadTransform('transformComposite.h5')
tp2TestRegImg = sitk.Resample(transImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp2TestReg = sitk.GetArrayFromImage(tp2TestRegImg)



In [35]:

    
visVolDiff(preproc(tp1Test), preproc(tp2TestReg))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0

Rigid Results

The registration definately moves us in the right direction; note the blobs on the bttom of the image (appx 900,500) are much closer than before. The transformation here, I think, is simply not strong enough to align the images. This could be due to a host of reasons:

Step size too small
Only 1 iteration of rigid alignment
Too few samples between images
Too much noise in input images

Next step is increasing the number of samples. Hopefully this will give the alg a better idea of the large structures of the image, and improve registration.

On the bright side, registration of 5 slices took about 30 seconds, wehereas it took SUBSTANTIALLY longer before.



In [41]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['20']
pMap['RequiredRatioOfValidSamples'] = ['.15']

SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 27.4 s, sys: 344 ms, total: 27.7 s
Wall time: 5.81 s



In [42]:

    
transformParameterMap = SimpleElastix.GetTransformParameterMap()

transImg = sitk.GetImageFromArray(tp2Test)

imageFilter = sitk.ReadTransform('transformComposite.h5')
tp2TestRegImg = sitk.Resample(transImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp2TestReg = sitk.GetArrayFromImage(tp2TestRegImg)

visVolDiff(preproc(tp1Test), preproc(tp2TestReg))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0



In [11]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['30']
pMap['RequiredRatioOfValidSamples'] = ['.35']

SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0






    



---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-11-54717cd49ab4> in <module>()
----> 1 get_ipython().run_cell_magic(u'time', u'', u"SimpleElastix = sitk.SimpleElastix()\n\nSimpleElastix.LogToConsoleOn()\n\nimg2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))\nnib.save(img2, 'fixed.nii')\nimg3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))\nnib.save(img3, 'moving.nii')\n\nSimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))\nSimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))\n\npMap = sitk.GetDefaultParameterMap('rigid')\npMap['AutomaticTransformInitialization'] =['true']\npMap['MaximumNumberOfSamplingAttempts'] = ['30']\npMap['RequiredRatioOfValidSamples'] = ['.35']\n\nSimpleElastix.SetParameterMap(pMap)\nSimpleElastix.Execute()")

/usr/local/lib/python2.7/dist-packages/IPython/core/interactiveshell.pyc in run_cell_magic(self, magic_name, line, cell)
   2113             magic_arg_s = self.var_expand(line, stack_depth)
   2114             with self.builtin_trap:
-> 2115                 result = fn(magic_arg_s, cell)
   2116             return result
   2117 

<decorator-gen-59> in time(self, line, cell, local_ns)

/usr/local/lib/python2.7/dist-packages/IPython/core/magic.pyc in <lambda>(f, *a, **k)
    186     # but it's overkill for just that one bit of state.
    187     def magic_deco(arg):
--> 188         call = lambda f, *a, **k: f(*a, **k)
    189 
    190         if callable(arg):

/usr/local/lib/python2.7/dist-packages/IPython/core/magics/execution.pyc in time(self, line, cell, local_ns)
   1178         else:
   1179             st = clock2()
-> 1180             exec(code, glob, local_ns)
   1181             end = clock2()
   1182             out = None

<timed exec> in <module>()

/usr/local/lib/python2.7/dist-packages/SimpleITK/SimpleITK.pyc in Execute(self)
  10028     def Execute(self):
  10029         """Execute(SimpleElastix self) -> Image"""
> 10030         return _SimpleITK.SimpleElastix_Execute(self)
  10031 
  10032 

RuntimeError: Exception thrown in SimpleITK SimpleElastix_Execute: /home/bstadt/simpleElastix/SimpleElastix/Code/Elastix/src/sitkSimpleElastixImpl.cxx:208:
sitk::ERROR: 
itk::ExceptionObject (0x59fd940)
Location: "unknown" 
File: /home/bstadt/simpleElastix/build/elastix/src/Core/Main/elxElastixFilter.hxx
Line: 240
Description: itk::ERROR: Self(0x5a65860): Internal elastix error: See elastix log (use LogToConsoleOn() or LogToFileOn()).

Looks like pumping up the number of samples only gets us so far before simple elastix can't handle it.



In [9]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('affine')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['20']
pMap['RequiredRatioOfValidSamples'] = ['.05']

SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 27.9 s, sys: 812 ms, total: 28.7 s
Wall time: 5.98 s



In [10]:

    
transformParameterMap = SimpleElastix.GetTransformParameterMap()

transImg = sitk.GetImageFromArray(tp2Test)

imageFilter = sitk.ReadTransform('transformComposite.h5')
tp2TestRegImg = sitk.Resample(transImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp2TestReg = sitk.GetArrayFromImage(tp2TestRegImg)

visVolDiff(preproc(tp1Test), preproc(tp2TestReg))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0

Affine gets us just about as close as rigid. This tells me that neither algorithm has enough iterations to converge - so time to figure out a way to increase that

Time is still looking good though



In [12]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['20']
pMap['RequiredRatioOfValidSamples'] = ['.05']
pMap['MaximumNumberOfIterations'] = ['100']
SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 16.8 s, sys: 348 ms, total: 17.2 s
Wall time: 4.03 s



In [13]:

    
transformParameterMap = SimpleElastix.GetTransformParameterMap()

transImg = sitk.GetImageFromArray(tp2Test)

imageFilter = sitk.ReadTransform('transformComposite.h5')
tp2TestRegImg = sitk.Resample(transImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp2TestReg = sitk.GetArrayFromImage(tp2TestRegImg)

visVolDiff(preproc(tp1Test), preproc(tp2TestReg))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0



In [14]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['20']
pMap['RequiredRatioOfValidSamples'] = ['.05']
pMap['MaximumNumberOfIterations'] = ['500']
SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 47 s, sys: 692 ms, total: 47.7 s
Wall time: 9.18 s



In [15]:

    
transformParameterMap = SimpleElastix.GetTransformParameterMap()

transImg = sitk.GetImageFromArray(tp2Test)

imageFilter = sitk.ReadTransform('transformComposite.h5')
tp2TestRegImg = sitk.Resample(transImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp2TestReg = sitk.GetArrayFromImage(tp2TestRegImg)

visVolDiff(preproc(tp1Test), preproc(tp2TestReg))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0

Number of iterations really doesnt change much - this is indicitive of the solver getting stuck in a local minimum. Im going to try to increase the step size in some way.



In [18]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['20']
pMap['RequiredRatioOfValidSamples'] = ['.05']
pMap['MaximumNumberOfIterations'] = ['100']
pMap['MaximumStepLength'] = ['5.0']
pMap['UseAdaptiveStepSizes'] = ['false']
SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0






    



---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-18-17d19d179731> in <module>()
----> 1 get_ipython().run_cell_magic(u'time', u'', u"SimpleElastix = sitk.SimpleElastix()\n\nSimpleElastix.LogToConsoleOn()\n\nimg2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))\nnib.save(img2, 'fixed.nii')\nimg3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))\nnib.save(img3, 'moving.nii')\n\nSimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))\nSimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))\n\npMap = sitk.GetDefaultParameterMap('rigid')\npMap['AutomaticTransformInitialization'] =['true']\npMap['MaximumNumberOfSamplingAttempts'] = ['20']\npMap['RequiredRatioOfValidSamples'] = ['.05']\npMap['MaximumNumberOfIterations'] = ['100']\npMap['MaximumStepLength'] = ['5.0']\npMap['UseAdaptiveStepSizes'] = ['false']\nSimpleElastix.SetParameterMap(pMap)\nSimpleElastix.Execute()")

/usr/local/lib/python2.7/dist-packages/IPython/core/interactiveshell.pyc in run_cell_magic(self, magic_name, line, cell)
   2113             magic_arg_s = self.var_expand(line, stack_depth)
   2114             with self.builtin_trap:
-> 2115                 result = fn(magic_arg_s, cell)
   2116             return result
   2117 

<decorator-gen-59> in time(self, line, cell, local_ns)

/usr/local/lib/python2.7/dist-packages/IPython/core/magic.pyc in <lambda>(f, *a, **k)
    186     # but it's overkill for just that one bit of state.
    187     def magic_deco(arg):
--> 188         call = lambda f, *a, **k: f(*a, **k)
    189 
    190         if callable(arg):

/usr/local/lib/python2.7/dist-packages/IPython/core/magics/execution.pyc in time(self, line, cell, local_ns)
   1178         else:
   1179             st = clock2()
-> 1180             exec(code, glob, local_ns)
   1181             end = clock2()
   1182             out = None

<timed exec> in <module>()

/usr/local/lib/python2.7/dist-packages/SimpleITK/SimpleITK.pyc in Execute(self)
  10028     def Execute(self):
  10029         """Execute(SimpleElastix self) -> Image"""
> 10030         return _SimpleITK.SimpleElastix_Execute(self)
  10031 
  10032 

RuntimeError: Exception thrown in SimpleITK SimpleElastix_Execute: /home/bstadt/simpleElastix/SimpleElastix/Code/Elastix/src/sitkSimpleElastixImpl.cxx:208:
sitk::ERROR: 
itk::ExceptionObject (0x5ba1390)
Location: "unknown" 
File: /home/bstadt/simpleElastix/build/elastix/src/Core/Main/elxElastixFilter.hxx
Line: 240
Description: itk::ERROR: Self(0xbb6da40): Internal elastix error: See elastix log (use LogToConsoleOn() or LogToFileOn()).

Changing the step length keeps crashing SimpleITK... I guess that means that It isn't a valid approach. Maybe, instead, change preprocessing so there is less noise for the registration to grab on to?



In [7]:

    
def preproc2(img):
    binImg = adaptiveThreshold(img, 5, 5)
    outImg = np.stack([cv2.erode(sub, None, 1) for sub in binImg])
    return outImg



In [16]:

    
visVolDiff(preproc2(tp1Test), preproc2(tp2Test))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0

Ok super much less noise for the algorthm to grab on to, with a bunch of clear landmarks still visible.



In [17]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc2(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc2(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['200']
pMap['RequiredRatioOfValidSamples'] = ['.05']
pMap['MaximumNumberOfIterations'] = ['100']
SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 14.4 s, sys: 620 ms, total: 15 s
Wall time: 3.69 s



In [19]:

    
transformParameterMap = SimpleElastix.GetTransformParameterMap()

transImg = sitk.GetImageFromArray(tp2Test)

imageFilter = sitk.ReadTransform('transformComposite.h5')
tp2TestRegImg = sitk.Resample(transImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp2TestReg = sitk.GetArrayFromImage(tp2TestRegImg)

visVolDiff(preproc2(tp1Test), preproc2(tp2TestReg))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0

Ok so I got super confused at this point for a while. The upper right hand corner of the image does give a clue though - notice how the green is moved to surround the red? Maybe this is the optimal solution for mutul information loss, but not for us. If I change the loss to l2, maybe things will turn out better. Worth a shot.



In [8]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['200']
pMap['CheckNumberOfSamples'] = ['false']
pMap['MaximumNumberOfIterations'] = ['100']
pMap['Metric'] = ['AdvancedNormalizedCorrelation']
SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 14 s, sys: 668 ms, total: 14.7 s
Wall time: 3.21 s



In [10]:

    
transformParameterMap = SimpleElastix.GetTransformParameterMap()

transImg = sitk.GetImageFromArray(tp2Test)

imageFilter = sitk.ReadTransform('transformComposite.h5')
tp2TestRegImg = sitk.Resample(transImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp2TestReg = sitk.GetArrayFromImage(tp2TestRegImg)

visVolDiff(preproc(tp1Test), preproc(tp2TestReg))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0



In [14]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['200']
pMap['CheckNumberOfSamples'] = ['false']
pMap['MaximumNumberOfIterations'] = ['100']
pMap['NumberOfSpatialSamples'] = ['1000000']
pMap['Metric'] = ['AdvancedNormalizedCorrelation']
SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 3min 2s, sys: 396 ms, total: 3min 2s
Wall time: 1min 15s



In [16]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc2(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc2(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['200']
pMap['RequiredRatioOfValidSamples'] = ['.05']
pMap['MaximumNumberOfIterations'] = ['500']

pMap['Metric'] = ['AdvancedMeanSquares']

SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 27.7 s, sys: 788 ms, total: 28.5 s
Wall time: 5.78 s



In [22]:

    
transformParameterMap = SimpleElastix.GetTransformParameterMap()

transImg = sitk.GetImageFromArray(tp2Test)
baseImg = sitk.GetImageFromArray(tp1Test)

imageFilter = sitk.ReadTransform('transformComposite.h5')
tp2TestRegImg = sitk.Resample(transImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp2TestReg = sitk.GetArrayFromImage(tp2TestRegImg)

imageFilter = sitk.ReadTransform('transformInverseComposite.h5')
tp1TestRegImg = sitk.Resample(baseImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp1TestReg = sitk.GetArrayFromImage(tp1TestRegImg)

visVolDiff(preproc2(tp1TestReg), preproc2(tp2TestReg))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0

Progress!

Looks like changing to l2 loss is definately the way to go. This, combined with the preprocessing erosion, has yielded the best results so far. And the entire transform remains fast at about a second per slice. Time to amp up the number of samples to see if precision increases



In [23]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc2(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc2(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['200']
pMap['RequiredRatioOfValidSamples'] = ['.05']
pMap['MaximumNumberOfIterations'] = ['500']
pMap['NumberOfSpatialSamples'] = ['100000']
pMap['Metric'] = ['AdvancedMeanSquares']

SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 1min 46s, sys: 1.52 s, total: 1min 48s
Wall time: 47.9 s



In [24]:

    
transformParameterMap = SimpleElastix.GetTransformParameterMap()

transImg = sitk.GetImageFromArray(tp2Test)
baseImg = sitk.GetImageFromArray(tp1Test)

imageFilter = sitk.ReadTransform('transformComposite.h5')
tp2TestRegImg = sitk.Resample(transImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp2TestReg = sitk.GetArrayFromImage(tp2TestRegImg)

imageFilter = sitk.ReadTransform('transformInverseComposite.h5')
tp1TestRegImg = sitk.Resample(baseImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp1TestReg = sitk.GetArrayFromImage(tp1TestRegImg)

visVolDiff(preproc2(tp1TestReg), preproc2(tp2TestReg))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0



In [25]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc2(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc2(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('affine')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['200']
pMap['RequiredRatioOfValidSamples'] = ['.05']
pMap['MaximumNumberOfIterations'] = ['500']
pMap['NumberOfSpatialSamples'] = ['100000']
pMap['Metric'] = ['AdvancedMeanSquares']

SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 1min 46s, sys: 1.61 s, total: 1min 48s
Wall time: 48.4 s



In [32]:

    
transformParameterMap = SimpleElastix.GetTransformParameterMap()

transImg = sitk.GetImageFromArray(tp2Test)
baseImg = sitk.GetImageFromArray(tp1Test)

imageFilter = sitk.ReadTransform('transformComposite.h5')
tp2TestRegImg = sitk.Resample(transImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp2TestReg = sitk.GetArrayFromImage(tp2TestRegImg)

imageFilter = sitk.ReadTransform('transformInverseComposite.h5')
tp1TestRegImg = sitk.Resample(baseImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp1TestReg = sitk.GetArrayFromImage(tp1TestRegImg)

visVolDiff(preproc2(tp1TestReg), preproc2(tp2TestReg))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0



In [55]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc2(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc2(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('bspline')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['500']
pMap['RequiredRatioOfValidSamples'] = ['.05']
pMap['MaximumNumberOfIterations'] = ['500']
pMap['Metric'] = ['AdvancedMeanSquares']

SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 2min 36s, sys: 1.8 s, total: 2min 37s
Wall time: 1min 38s



In [57]:

    
transformParameterMap = SimpleElastix.GetTransformParameterMap()

transImg = sitk.GetImageFromArray(tp2Test)
baseImg = sitk.GetImageFromArray(tp1Test)

imageFilter = sitk.ReadTransform('transformComposite.h5')
tp2TestRegImg = sitk.Resample(transImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp2TestReg = sitk.GetArrayFromImage(tp2TestRegImg)

imageFilter = sitk.ReadTransform('transformInverseComposite.h5')
tp1TestRegImg = sitk.Resample(baseImg,imageFilter, sitk.sitkGaussian, sitk.sitkFloat32)
tp1TestReg = sitk.GetArrayFromImage(tp1TestRegImg)

visVolDiff(preproc2(tp1TestReg), preproc2(tp2TestReg))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0



In [58]:

    
outImg = sitk.GetArrayFromImage(SimpleElastix.GetResultImage())



In [59]:

    
out = np.rollaxis(outImg, 2, 0)



In [60]:

    
plt.figure()
plt.imshow(out[2], cmap='gray')
plt.show()

This is awkward

It seems like through the testing I have not been extracting the transformations correctly. Above is the reported disperity from the transform, where white indicates unregistered portions of the image.

See how little white there is? And see how my version has a bunch of disperity? That means I am extracting the transformation wrong.

On the bright side, this means that things are working!

Time to retry rigid with the same stratgy



In [17]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['200']
pMap['RequiredRatioOfValidSamples'] = ['.05']
pMap['WriteCompositeTransform'] = ['true']

SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 27.3 s, sys: 760 ms, total: 28 s
Wall time: 5.8 s



In [62]:

    
outImg = sitk.GetArrayFromImage(SimpleElastix.GetResultImage())
out = np.rollaxis(outImg, 2, 0)
plt.figure()
plt.imshow(out[2], cmap='gray')
plt.show()



In [25]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc(tp1Test), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc(tp2Test), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['200']
pMap['RequiredRatioOfValidSamples'] = ['.05']
pMap['WriteCompositeTransform'] = ['true']

SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
CPU times: user 28 s, sys: 660 ms, total: 28.6 s
Wall time: 6 s



In [41]:

    
transformParameterMap = SimpleElastix.GetTransformParameterMap()

imgFilter = sitk.SimpleTransformix()
imgFilter.SetTransformParameterMap(transformParameterMap)
imgFilter.SetMovingImage(sitk.GetImageFromArray(tp1Test))
imgFilter.Execute()

tp1RegImg = imgFilter.GetResultImage()
tp1TestReg = sitk.GetArrayFromImage(tp1RegImg)
tp1TestReg = np.rollaxis(tp1TestReg, 2, 0)
visVolDiff(preproc(tp1TestReg), preproc(tp2Test))









    



0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0



In [43]:

    
%%time 
SimpleElastix = sitk.SimpleElastix()

SimpleElastix.LogToConsoleOn()

img2 = nib.Nifti1Image(preproc(tp1ChanList[1]), np.eye(4))
nib.save(img2, 'fixed.nii')
img3 = nib.Nifti1Image(preproc(tp2ChanList[1]), np.eye(4))
nib.save(img3, 'moving.nii')

SimpleElastix.SetFixedImage(sitk.ReadImage('fixed.nii'))
SimpleElastix.SetMovingImage(sitk.ReadImage('moving.nii'))

pMap = sitk.GetDefaultParameterMap('rigid')
pMap['AutomaticTransformInitialization'] =['true']
pMap['MaximumNumberOfSamplingAttempts'] = ['200']
pMap['RequiredRatioOfValidSamples'] = ['.05']
pMap['WriteCompositeTransform'] = ['true']

SimpleElastix.SetParameterMap(pMap)
SimpleElastix.Execute()









    



0.00357142857143
0.00714285714286
0.0107142857143
0.0142857142857
0.0178571428571
0.0214285714286
0.025
0.0285714285714
0.0321428571429
0.0357142857143
0.0392857142857
0.0428571428571
0.0464285714286
0.05
0.0535714285714
0.0571428571429
0.0607142857143
0.0642857142857
0.0678571428571
0.0714285714286
0.075
0.0785714285714
0.0821428571429
0.0857142857143
0.0892857142857
0.0928571428571
0.0964285714286
0.1
0.103571428571
0.107142857143
0.110714285714
0.114285714286
0.117857142857
0.121428571429
0.125
0.128571428571
0.132142857143
0.135714285714
0.139285714286
0.142857142857
0.146428571429
0.15
0.153571428571
0.157142857143
0.160714285714
0.164285714286
0.167857142857
0.171428571429
0.175
0.178571428571
0.182142857143
0.185714285714
0.189285714286
0.192857142857
0.196428571429
0.2
0.203571428571
0.207142857143
0.210714285714
0.214285714286
0.217857142857
0.221428571429
0.225
0.228571428571
0.232142857143
0.235714285714
0.239285714286
0.242857142857
0.246428571429
0.25
0.253571428571
0.257142857143
0.260714285714
0.264285714286
0.267857142857
0.271428571429
0.275
0.278571428571
0.282142857143
0.285714285714
0.289285714286
0.292857142857
0.296428571429
0.3
0.303571428571
0.307142857143
0.310714285714
0.314285714286
0.317857142857
0.321428571429
0.325
0.328571428571
0.332142857143
0.335714285714
0.339285714286
0.342857142857
0.346428571429
0.35
0.353571428571
0.357142857143
0.360714285714
0.364285714286
0.367857142857
0.371428571429
0.375
0.378571428571
0.382142857143
0.385714285714
0.389285714286
0.392857142857
0.396428571429
0.4
0.403571428571
0.407142857143
0.410714285714
0.414285714286
0.417857142857
0.421428571429
0.425
0.428571428571
0.432142857143
0.435714285714
0.439285714286
0.442857142857
0.446428571429
0.45
0.453571428571
0.457142857143
0.460714285714
0.464285714286
0.467857142857
0.471428571429
0.475
0.478571428571
0.482142857143
0.485714285714
0.489285714286
0.492857142857
0.496428571429
0.5
0.503571428571
0.507142857143
0.510714285714
0.514285714286
0.517857142857
0.521428571429
0.525
0.528571428571
0.532142857143
0.535714285714
0.539285714286
0.542857142857
0.546428571429
0.55
0.553571428571
0.557142857143
0.560714285714
0.564285714286
0.567857142857
0.571428571429
0.575
0.578571428571
0.582142857143
0.585714285714
0.589285714286
0.592857142857
0.596428571429
0.6
0.603571428571
0.607142857143
0.610714285714
0.614285714286
0.617857142857
0.621428571429
0.625
0.628571428571
0.632142857143
0.635714285714
0.639285714286
0.642857142857
0.646428571429
0.65
0.653571428571
0.657142857143
0.660714285714
0.664285714286
0.667857142857
0.671428571429
0.675
0.678571428571
0.682142857143
0.685714285714
0.689285714286
0.692857142857
0.696428571429
0.7
0.703571428571
0.707142857143
0.710714285714
0.714285714286
0.717857142857
0.721428571429
0.725
0.728571428571
0.732142857143
0.735714285714
0.739285714286
0.742857142857
0.746428571429
0.75
0.753571428571
0.757142857143
0.760714285714
0.764285714286
0.767857142857
0.771428571429
0.775
0.778571428571
0.782142857143
0.785714285714
0.789285714286
0.792857142857
0.796428571429
0.8
0.803571428571
0.807142857143
0.810714285714
0.814285714286
0.817857142857
0.821428571429
0.825
0.828571428571
0.832142857143
0.835714285714
0.839285714286
0.842857142857
0.846428571429
0.85
0.853571428571
0.857142857143
0.860714285714
0.864285714286
0.867857142857
0.871428571429
0.875
0.878571428571
0.882142857143
0.885714285714
0.889285714286
0.892857142857
0.896428571429
0.9
0.903571428571
0.907142857143
0.910714285714
0.914285714286
0.917857142857
0.921428571429
0.925
0.928571428571
0.932142857143
0.935714285714
0.939285714286
0.942857142857
0.946428571429
0.95
0.953571428571
0.957142857143
0.960714285714
0.964285714286
0.967857142857
0.971428571429
0.975
0.978571428571
0.982142857143
0.985714285714
0.989285714286
0.992857142857
0.996428571429
1.0
0.00357142857143
0.00714285714286
0.0107142857143
0.0142857142857
0.0178571428571
0.0214285714286
0.025
0.0285714285714
0.0321428571429
0.0357142857143
0.0392857142857
0.0428571428571
0.0464285714286
0.05
0.0535714285714
0.0571428571429
0.0607142857143
0.0642857142857
0.0678571428571
0.0714285714286
0.075
0.0785714285714
0.0821428571429
0.0857142857143
0.0892857142857
0.0928571428571
0.0964285714286
0.1
0.103571428571
0.107142857143
0.110714285714
0.114285714286
0.117857142857
0.121428571429
0.125
0.128571428571
0.132142857143
0.135714285714
0.139285714286
0.142857142857
0.146428571429
0.15
0.153571428571
0.157142857143
0.160714285714
0.164285714286
0.167857142857
0.171428571429
0.175
0.178571428571
0.182142857143
0.185714285714
0.189285714286
0.192857142857
0.196428571429
0.2
0.203571428571
0.207142857143
0.210714285714
0.214285714286
0.217857142857
0.221428571429
0.225
0.228571428571
0.232142857143
0.235714285714
0.239285714286
0.242857142857
0.246428571429
0.25
0.253571428571
0.257142857143
0.260714285714
0.264285714286
0.267857142857
0.271428571429
0.275
0.278571428571
0.282142857143
0.285714285714
0.289285714286
0.292857142857
0.296428571429
0.3
0.303571428571
0.307142857143
0.310714285714
0.314285714286
0.317857142857
0.321428571429
0.325
0.328571428571
0.332142857143
0.335714285714
0.339285714286
0.342857142857
0.346428571429
0.35
0.353571428571
0.357142857143
0.360714285714
0.364285714286
0.367857142857
0.371428571429
0.375
0.378571428571
0.382142857143
0.385714285714
0.389285714286
0.392857142857
0.396428571429
0.4
0.403571428571
0.407142857143
0.410714285714
0.414285714286
0.417857142857
0.421428571429
0.425
0.428571428571
0.432142857143
0.435714285714
0.439285714286
0.442857142857
0.446428571429
0.45
0.453571428571
0.457142857143
0.460714285714
0.464285714286
0.467857142857
0.471428571429
0.475
0.478571428571
0.482142857143
0.485714285714
0.489285714286
0.492857142857
0.496428571429
0.5
0.503571428571
0.507142857143
0.510714285714
0.514285714286
0.517857142857
0.521428571429
0.525
0.528571428571
0.532142857143
0.535714285714
0.539285714286
0.542857142857
0.546428571429
0.55
0.553571428571
0.557142857143
0.560714285714
0.564285714286
0.567857142857
0.571428571429
0.575
0.578571428571
0.582142857143
0.585714285714
0.589285714286
0.592857142857
0.596428571429
0.6
0.603571428571
0.607142857143
0.610714285714
0.614285714286
0.617857142857
0.621428571429
0.625
0.628571428571
0.632142857143
0.635714285714
0.639285714286
0.642857142857
0.646428571429
0.65
0.653571428571
0.657142857143
0.660714285714
0.664285714286
0.667857142857
0.671428571429
0.675
0.678571428571
0.682142857143
0.685714285714
0.689285714286
0.692857142857
0.696428571429
0.7
0.703571428571
0.707142857143
0.710714285714
0.714285714286
0.717857142857
0.721428571429
0.725
0.728571428571
0.732142857143
0.735714285714
0.739285714286
0.742857142857
0.746428571429
0.75
0.753571428571
0.757142857143
0.760714285714
0.764285714286
0.767857142857
0.771428571429
0.775
0.778571428571
0.782142857143
0.785714285714
0.789285714286
0.792857142857
0.796428571429
0.8
0.803571428571
0.807142857143
0.810714285714
0.814285714286
0.817857142857
0.821428571429
0.825
0.828571428571
0.832142857143
0.835714285714
0.839285714286
0.842857142857
0.846428571429
0.85
0.853571428571
0.857142857143
0.860714285714
0.864285714286
0.867857142857
0.871428571429
0.875
0.878571428571
0.882142857143
0.885714285714
0.889285714286
0.892857142857
0.896428571429
0.9
0.903571428571
0.907142857143
0.910714285714
0.914285714286
0.917857142857
0.921428571429
0.925
0.928571428571
0.932142857143
0.935714285714
0.939285714286
0.942857142857
0.946428571429
0.95
0.953571428571
0.957142857143
0.960714285714
0.964285714286
0.967857142857
0.971428571429
0.975
0.978571428571
0.982142857143
0.985714285714
0.989285714286
0.992857142857
0.996428571429
1.0
CPU times: user 12min 41s, sys: 49.6 s, total: 13min 30s
Wall time: 3min 4s



In [46]:

    
outImg = sitk.GetArrayFromImage(SimpleElastix.GetResultImage())
out = np.rollaxis(outImg, 2, 0)
for idx in range(0, out.shape[0], 25):
    plt.figure()
    plt.imshow(out[idx], cmap='gray')
    plt.title('Disperity')
    plt.show()

    plt.figure()
    plt.imshow(tp1ChanList[1][idx], cmap='gray')
    plt.title('Tp1')
    plt.show()
    
    plt.figure()
    plt.imshow(tp2ChanList[1][idx], cmap='gray')
    plt.title('Tp2')
    plt.show()
    
    print '\n'
    print '\n'

    print '\n'









    





<matplotlib.figure.Figure at 0x7ff941eca2d0>



In [ ]: