Kwame Comment Tests

I met with Kwame again on 1/23, and he was able to come up with a few suggestions that may help improve registration accuracy. His comments are briefly explored here



In [1]:

    
import sys
sys.path.insert(0, '../../../ndreg/')
sys.path.insert(0,'../code/functions/')

import ndreg
import math
import cv2
import pickle

import numpy as np
import SimpleITK as itk
import matplotlib.pyplot as plt
import plosLib as pLib
import connectLib as cLib
import connectLib as cLib

from affine import Affine
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
from random import randrange as rand
from random import uniform as floatRand
from PIL import Image



In [19]:

    
def register(inImg, baseImg, iterations, MI=False):
    return ndreg.imgAffineComposite(inImg,
                           baseImg,
                           iterations=iterations,
                           useMI=MI,
                           verbose=False)



In [3]:

    
def getTransform(pitchRange = .1,
                 rollRange = .1,
                 yawRange = .1,
                 xRange = 0.,
                 yRange = 0.,
                 zRange = 0.,
                 pitch=None,
                 yaw=None,
                 roll=None,
                 zT=None,
                 yT=None,
                 xT=None):
    
    #generate a random rigid body transform in given bounds
    a = floatRand(pitchRange-(pitchRange/2.), pitchRange+(pitchRange/2.))
    b = floatRand(yawRange-(yawRange/2.), yawRange+(yawRange/2.))
    c = floatRand(rollRange-(rollRange/2.), rollRange+(rollRange/2.))

    xt = floatRand(xRange-(xRange/2.), xRange+(xRange/2.))
    yt = floatRand(yRange-(yRange/2.), yRange+(yRange/2.))
    zt = floatRand(zRange-(zRange/2.), zRange+(zRange/2.))

    #set default params, if passed
    if yaw is not None:
        a = yaw
    if pitch is not None:
        b = pitch
    if roll is not None:
        c = roll 
    if xT is not None:
        xt = xT
    if yT is not None:
        yt = yT
    if zT is not None:
        zt = zT

    #generate the transform
    transform = np.stack([
        [math.cos(a)*math.cos(b), math.cos(a)*math.sin(b)*math.sin(c)-math.sin(a)*math.cos(c), math.cos(a)*math.sin(b)*math.cos(c)+math.sin(a)*math.sin(c), xt],
        [math.sin(a)*math.cos(b), math.sin(a)*math.sin(b)*math.sin(c)+math.cos(a)*math.cos(c), math.sin(a)*math.sin(b)*math.cos(c)-math.cos(a)*math.sin(c), yt],
        [-math.sin(b), math.cos(b)*math.sin(c), math.cos(b)*math.cos(c), zt],
        [0., 0., 0., 1.]
    ])
    return transform
    
    
def applyRigid(initialVolume, transform, verbose=False):
    #create a body to hold return values
    rigidMatrix = np.zeros_like(initialVolume)  
    
    #create vars to keep track of out of bounds percent
    out = 0.
    
    #convert every voxel
    for z in range(initialVolume.shape[0]):
        for y in range(initialVolume.shape[1]):
            for x in range(initialVolume.shape[2]):
                
                #shift the volume so that the transform is applied to the center, and not the top left back edge
                shift = [initialVolume.shape[2]/2., initialVolume.shape[1]/2., initialVolume.shape[0]/2.]
                shiftedPoint = list(np.subtract([x, y, z], shift))
                shiftedPoint.append(1.)
                
                #calculate the transform, drop the irrelavent 1 at the end
                new = np.dot(transform, shiftedPoint)[:3]
                
                #shift the volume back to its original spot after transform
                final = np.add(new, shift)
                
                #attempt to assign new spot
                try:
                    rigidMatrix[int(final[2]), int(final[1]), int(final[0])] = initialVolume[z, y, x]
                    
                #if transformed place is out of bounds, dont assign it
                except IndexError:
                        out+=1.
                        continue
    #print information, if required
    if verbose:
        print 'Out of bounds fraction: ', out/float(initialVolume.shape[0]*initialVolume.shape[1]*initialVolume.shape[2])
    
    return rigidMatrix



In [4]:

    
realData = pickle.load(open('../code/tests/synthDat/realDataRaw_t0.synth', 'r'))
print realData.shape









    



(280, 1024, 1024)



In [5]:

    
testData = realData[:7]
print testData.shape









    



(7, 1024, 1024)

Spacing Test

The goal of the spacing test is the determine if there is an effect of setting image spacing on NDReg outputs



In [10]:

    
def executeTest(testImg, baseImg, gradientIterations = 200, verbose=True):
    
    initialErr = np.sum((np.subtract(baseImg, testImg))**2)
    
    testImgITK = itk.GetImageFromArray(testImg)
    testImgITK.SetSpacing([.12, .12, .5])
    baseImgITK = itk.GetImageFromArray(baseImg)
    baseImgITK.SetSpacing([.12, .12, .5])
    
    if verbose:
        print '\tInitial Error: ', initialErr
    transform = register(testImgITK, baseImgITK, gradientIterations)
    
    if verbose:
        print '\tTransform:'
        print transform
         
    regImg = itk.GetArrayFromImage(
        ndreg.imgApplyAffine(
            testImgITK,
            transform,
            size=testImgITK.GetSize(),
            spacing=testImgITK.GetSpacing()
            )
        )

    finalErr = np.sum((np.subtract(baseImg, regImg))**2)
    errReductionRatio = np.subtract(initialErr, finalErr)/float(initialErr)
    
    if verbose:
        print '\tFinal Error: ', finalErr
        print '\tError Reduction Ratio: ', errReductionRatio
    
    return regImg, initialErr, finalErr, errReductionRatio



In [7]:

    
outVolume = cLib.otsuVox(pLib.pipeline(testData))



In [8]:

    
#As expected, slices 0 and 6 are ill defined, so remove them
outVolume = outVolume[1:6]
print outVolume.shape









    



(5, 1024, 1024)



In [9]:

    
#generate and apply a basic transform to the output
transform = getTransform(pitch=0., yaw=.15, roll=0., xT=0., yT=0., zT=0.)
transformVolume = applyRigid(outVolume, transform, True)









    



Out of bounds fraction:  0.0318670272827



In [20]:

    
#Check if NDReg handles the basic case
regImg, initialErr, finalErr, reductionRatio = executeTest(transformVolume, outVolume, 200, verbose=True)









    



	Initial Error:  110266.0
	Transform:
[0.9975932877997035, 0.05842507161859518, -0.05792548260651925, -0.06162104745726014, 0.9975625945148371, -0.055072088456448946, 0.054513273971632206, 0.05845168451847105, 0.9977840056493715, 0.04733375669943689, -0.053518760362313074, 0.07005213287769078]
	Final Error:  60874.0
	Error Reduction Ratio:  0.447934993561

Spacing Test Results

Error reduction ratio and Final Error are the same in this test as the baseline

Affine Test

The goal of this test is to see if changing the registration method from a composite function to an affine only function will improve NDReg's capability to work on our current data



In [21]:

    
def executeTest(testImg, baseImg, gradientIterations = 200, verbose=True):
    
    initialErr = np.sum((np.subtract(baseImg, testImg))**2)
    
    testImgITK = itk.GetImageFromArray(testImg)
    baseImgITK = itk.GetImageFromArray(baseImg)
    
    if verbose:
        print '\tInitial Error: ', initialErr
    transform = register(testImgITK, baseImgITK, gradientIterations)
    
    if verbose:
        print '\tTransform:'
        print transform
         
    regImg = itk.GetArrayFromImage(
        ndreg.imgApplyAffine(
            testImgITK,
            transform,
            size=testImgITK.GetSize(),
            spacing=testImgITK.GetSpacing()
            )
        )

    finalErr = np.sum((np.subtract(baseImg, regImg))**2)
    errReductionRatio = np.subtract(initialErr, finalErr)/float(initialErr)
    
    if verbose:
        print '\tFinal Error: ', finalErr
        print '\tError Reduction Ratio: ', errReductionRatio
    
    return regImg, initialErr, finalErr, errReductionRatio



In [22]:

    
def register(inImg, baseImg, iterations, MI=False):
    return ndreg.imgAffine(inImg,
                           baseImg,
                           iterations=iterations,
                           useMI=MI,
                           verbose=False)



In [23]:

    
regImg, initialErr, finalErr, reductionRatio = executeTest(transformVolume, outVolume, 200, verbose=True)









    



	Initial Error:  110266.0
	Transform:
[1.0314913643177692, 0.023057167756149584, 0.00012169425691472045, -0.020255069538658053, 0.9710260065164376, -0.0001237295433937009, 0.06257475638247516, 0.05753484988843319, 0.9990004111255613, 5.280293545418143e-05, -5.970241331080096e-05, 7.817365375544298e-05]
	Final Error:  60874.0
	Error Reduction Ratio:  0.447934993561

Affine Function Test Conclusion

Results have not changed



In [ ]: