[Scipy-svn] r3929 - trunk/scipy/ndimage
scipy-svn@scip...
scipy-svn@scip...
Tue Feb 12 18:12:01 CST 2008
Author: chris.burns
Date: 2008-02-12 18:11:57 -0600 (Tue, 12 Feb 2008)
New Revision: 3929
Added:
trunk/scipy/ndimage/_registration.py
Removed:
trunk/scipy/ndimage/registration.py
Log:
Make registration a private module. Issue development warning on import.
Copied: trunk/scipy/ndimage/_registration.py (from rev 3928, trunk/scipy/ndimage/registration.py)
===================================================================
--- trunk/scipy/ndimage/registration.py 2008-02-13 00:00:26 UTC (rev 3928)
+++ trunk/scipy/ndimage/_registration.py 2008-02-13 00:11:57 UTC (rev 3929)
@@ -0,0 +1,457 @@
+import math
+import os
+import numpy as N
+import scipy.ndimage._register as R
+import scipy.special as SP
+import scipy.ndimage as NDI
+import scipy.optimize as OPT
+import time
+
+# Issue warning regarding heavy development status of this module
+import warnings
+_msg = "The registration code is under heavy development and therefore the \
+public API will change in the future. The NIPY group is actively working on \
+this code, and has every intention of generalizing this for the Scipy \
+community. Use this module minimally, if at all, until it this warning is \
+removed."
+warnings.warn(_msg, UserWarning)
+
+# TODO: Add docstrings for public functions in extension code.
+# Add docstrings to extension code.
+#from numpy.lib import add_newdoc
+#add_newdoc('scipy.ndimage._register', 'register_histogram',
+# """A joint histogram used for registration module.
+# """)
+
+
+# anatomical MRI to test with
+inputname = 'ANAT1_V0001.img'
+filename = os.path.join(os.path.split(__file__)[0], inputname)
+
+def remap_image(image, parm_vector):
+ M_inverse = get_inverse_mappings(parm_vector)
+ # allocate the zero image
+ remaped_image = load_blank_image()
+ imdata = build_structs(step=1)
+ # trilinear interpolation mapping. to be replaced with splines
+ R.register_linear_resample(image['data'], remaped_image['data'], M_inverse, imdata['step'])
+ return remaped_image
+
+def get_inverse_mappings(parm_vector):
+ # get the inverse mapping to rotate the G matrix to F space following registration
+ imdata = build_structs(step=1)
+ # inverse angles and translations
+ imdata['parms'][0] = -parm_vector[0]
+ imdata['parms'][1] = -parm_vector[1]
+ imdata['parms'][2] = -parm_vector[2]
+ imdata['parms'][3] = -parm_vector[3]
+ imdata['parms'][4] = -parm_vector[4]
+ imdata['parms'][5] = -parm_vector[5]
+ M_inverse = build_rotate_matrix(imdata['parms'])
+ return M_inverse
+
+def python_coreg(image1, image2, imdata, ftype=1, smimage=0, lite=0, smhist=0,
+ method='nmi', opt_method='powell'):
+ # image1 is imageF and image2 is imageG in SPM lingo
+ # get these from get_test_images for the debug work
+ start = time.time()
+ # smooth of the images
+ if smimage:
+ image_F_xyz1 = filter_image_3D(image1['data'], image1['fwhm'], ftype)
+ image1['data'] = image_F_xyz1
+ image_F_xyz2 = filter_image_3D(image2['data'], image2['fwhm'], ftype)
+ image2['data'] = image_F_xyz2
+ parm_vector = multires_registration(image1, image2, imdata, lite, smhist, method, opt_method)
+ stop = time.time()
+ print 'Total Optimizer Time is ', (stop-start)
+ return parm_vector
+
+def get_test_images(alpha=0.0, beta=0.0, gamma=0.0):
+ image1 = load_image()
+ image2 = load_blank_image()
+ imdata = build_structs(step=1)
+ # allow the G image to be rotated for testing
+ imdata['parms'][0] = alpha
+ imdata['parms'][1] = beta
+ imdata['parms'][2] = gamma
+ image1['fwhm'] = build_fwhm(image1['mat'], imdata['step'])
+ image2['fwhm'] = build_fwhm(image2['mat'], imdata['step'])
+ M = build_rotate_matrix(imdata['parms'])
+ R.register_linear_resample(image1['data'], image2['data'], M, imdata['step'])
+ return image1, image2, imdata
+
+def multires_registration(image1, image2, imdata, lite, smhist, method, opt_method):
+ ret_histo=0
+ # zero out the start parameter; but this may be set to large values
+ # if the head is out of range and well off the optimal alignment skirt
+ imdata['parms'][0:5] = 0.0
+ # make the step a scalar to can put in a multi-res loop
+ loop = range(imdata['sample'].size)
+ x = imdata['parms']
+ for i in loop:
+ step = imdata['sample'][i]
+ imdata['step'][:] = step
+ optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
+ p_args = (optfunc_args,)
+ if opt_method=='powell':
+ print 'POWELL multi-res registration step size ', step
+ print 'vector ', x
+ x = OPT.fmin_powell(optimize_function, x, args=p_args, callback=callback_powell)
+ elif opt_method=='cg':
+ print 'CG multi-res registration step size ', step
+ print 'vector ', x
+ x = OPT.fmin_cg(optimize_function, x, args=p_args, callback=callback_cg)
+ elif opt_method=='hybrid':
+ if i==0:
+ print 'Hybrid POWELL multi-res registration step size ', step
+ print 'vector ', x
+ lite = 0
+ optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
+ p_args = (optfunc_args,)
+ x = OPT.fmin_powell(optimize_function, x, args=p_args, callback=callback_powell)
+ elif i==1:
+ print 'Hybrid CG multi-res registration step size ', step
+ print 'vector ', x
+ lite = 1
+ optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
+ p_args = (optfunc_args,)
+ x = OPT.fmin_cg(optimize_function, x, args=p_args, callback=callback_cg)
+
+ return x
+
+
+def test_image_filter(image, imdata, ftype=2):
+ # test the 3D image filter on an image. ftype 1 is SPM, ftype 2 is simple Gaussian
+ image['fwhm'] = build_fwhm(image['mat'], imdata['step'])
+ filt_image = filter_image_3D(image['data'], image['fwhm'], ftype)
+ return filt_image
+
+def callback_powell(x):
+ print 'Parameter Vector from Powell: - '
+ print x
+ return
+
+def callback_cg(x):
+ print 'Parameter Vector from Conjugate Gradient: - '
+ print x
+ return
+
+def test_alignment(image1, image2, imdata, method='ncc', lite=0, smhist=0,
+ alpha=0.0, beta=0.0, gamma=0.0, ret_histo=0):
+
+ # to test the cost function and view the joint histogram
+ # for 2 images. used for debug
+ imdata['parms'][0] = alpha
+ imdata['parms'][1] = beta
+ imdata['parms'][2] = gamma
+ M = build_rotate_matrix(imdata['parms'])
+ optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
+
+ if ret_histo:
+ cost, joint_histogram = optimize_function(imdata['parms'], optfunc_args)
+ return cost, joint_histogram
+ else:
+ cost = optimize_function(imdata['parms'], optfunc_args)
+ return cost
+
+
+def smooth_kernel(fwhm, x, ktype=1):
+ eps = 0.00001
+ s = N.square((fwhm/math.sqrt(8.0*math.log(2.0)))) + eps
+ if ktype==1:
+ # from SPM: Gauss kernel convolved with 1st degree B spline
+ w1 = 0.5 * math.sqrt(2.0/s)
+ w2 = -0.5 / s
+ w3 = math.sqrt((s*math.pi) /2.0)
+ kernel = 0.5*(SP.erf(w1*(x+1))*(x+1) + SP.erf(w1*(x-1))*(x-1) - 2.0*SP.erf(w1*(x))*(x) +
+ w3*(N.exp(w2*N.square(x+1))) + N.exp(w2*(N.square(x-1))) - 2.0*N.exp(w2*N.square(x)))
+ kernel[kernel<0] = 0
+ kernel = kernel / kernel.sum()
+ else:
+ # Gauss kernel
+ kernel = (1.0/math.sqrt(2.0*math.pi*s)) * N.exp(-N.square(x)/(2.0*s))
+ kernel = kernel / kernel.sum()
+
+ return kernel
+
+def filter_image_3D(imageRaw, fwhm, ftype=2):
+ p = N.ceil(2*fwhm[0]).astype(int)
+ x = N.array(range(-p, p+1))
+ kernel_x = smooth_kernel(fwhm[0], x, ktype=ftype)
+ p = N.ceil(2*fwhm[1]).astype(int)
+ x = N.array(range(-p, p+1))
+ kernel_y = smooth_kernel(fwhm[1], x, ktype=ftype)
+ p = N.ceil(2*fwhm[2]).astype(int)
+ x = N.array(range(-p, p+1))
+ kernel_z = smooth_kernel(fwhm[2], x, ktype=ftype)
+ output=None
+ # 3D filter in 3 1D separable stages
+ axis = 0
+ image_F_x = NDI.correlate1d(imageRaw, kernel_x, axis, output)
+ axis = 1
+ image_F_xy = NDI.correlate1d(image_F_x, kernel_y, axis, output)
+ axis = 2
+ image_F_xyz = NDI.correlate1d(image_F_xy, kernel_z, axis, output)
+ return image_F_xyz
+
+
+def resample_image(smimage=0, ftype=2, alpha=0.0, beta=0.0, gamma=0.0,
+ Tx=0.0, Ty=0.0, Tz=0.0, stepsize=1):
+
+ # takes an image and 3D rotate using trilinear interpolation
+ image1 = load_image()
+ image2 = load_blank_image()
+ imdata = build_structs(step=stepsize)
+ imdata['parms'][0] = alpha
+ imdata['parms'][1] = beta
+ imdata['parms'][2] = gamma
+ imdata['parms'][3] = Tx
+ imdata['parms'][4] = Ty
+ imdata['parms'][5] = Tz
+ image1['fwhm'] = build_fwhm(image1['mat'], imdata['step'])
+ image2['fwhm'] = build_fwhm(image2['mat'], imdata['step'])
+ M = build_rotate_matrix(imdata['parms'])
+ if smimage:
+ image_F_xyz1 = filter_image_3D(image1['data'], image1['fwhm'], ftype)
+ image1['data'] = image_F_xyz1
+ image_F_xyz2 = filter_image_3D(image2['data'], image2['fwhm'], ftype)
+ image2['data'] = image_F_xyz2
+
+ R.register_linear_resample(image1['data'], image2['data'], M, imdata['step'])
+
+ return image2
+
+
+def build_fwhm(M, S):
+ view_3x3 = N.square(M[0:3, 0:3])
+ vxg = N.sqrt(view_3x3.sum(axis=0))
+ # assumes that sampling is the same for xyz
+ size = N.array([1,1,1])*S[0]
+ x = N.square(size) - N.square(vxg)
+ # clip
+ x[x<0] = 0
+ fwhm = N.sqrt(x) / vxg
+ # pathology when stepsize = 1 for MAT equal to the identity matrix
+ fwhm[fwhm==0] = 1
+ # return the 3D Gaussian kernel width (xyz)
+ return fwhm
+
+def load_image(imagename=filename, rows=256, cols=256, layers=90):
+ ImageVolume = N.fromfile(imagename, dtype=N.uint16).reshape(layers, rows, cols);
+ # clip to 8 bits. this will be rescale to 8 bits for fMRI
+ ImageVolume[ImageVolume>255] = 255
+ # voxel to pixel is identity for this simulation using anatomical MRI volume
+ # 4x4 matrix
+ M = N.eye(4, dtype=N.float64);
+ # dimensions
+ D = N.zeros(3, dtype=N.int32);
+ # Gaussian kernel - fill in with build_fwhm()
+ F = N.zeros(3, dtype=N.float64);
+ D[0] = rows
+ D[1] = cols
+ D[2] = layers
+ # make sure the data type is uchar
+ ImageVolume = ImageVolume.astype(N.uint8)
+ image = {'data' : ImageVolume, 'mat' : M, 'dim' : D, 'fwhm' : F}
+ return image
+
+def load_blank_image(rows=256, cols=256, layers=90):
+ ImageVolume = N.zeros(layers*rows*cols, dtype=N.uint16).reshape(layers, rows, cols);
+ # voxel to pixel is identity for this simulation using anatomical MRI volume
+ # 4x4 matrix
+ M = N.eye(4, dtype=N.float64);
+ # dimensions
+ D = N.zeros(3, dtype=N.int32);
+ # Gaussian kernel - fill in with build_fwhm()
+ F = N.zeros(3, dtype=N.float64);
+ D[0] = rows
+ D[1] = cols
+ D[2] = layers
+ # make sure the data type is uchar
+ ImageVolume = ImageVolume.astype(N.uint8)
+ image = {'data' : ImageVolume, 'mat' : M, 'dim' : D, 'fwhm' : F}
+ return image
+
+def optimize_function(x, optfunc_args):
+ image_F = optfunc_args[0]
+ image_G = optfunc_args[1]
+ sample_vector = optfunc_args[2]
+ fwhm = optfunc_args[3]
+ do_lite = optfunc_args[4]
+ smooth = optfunc_args[5]
+ method = optfunc_args[6]
+ ret_histo = optfunc_args[7]
+
+ rot_matrix = build_rotate_matrix(x)
+ cost = 0.0
+ epsilon = 2.2e-16
+ # image_G is base image
+ # image_F is the to-be-rotated image
+ # rot_matrix is the 4x4 constructed (current angles and translates) transform matrix
+ # sample_vector is the subsample vector for x-y-z
+
+ F_inv = N.linalg.inv(image_F['mat'])
+ composite = N.dot(F_inv, rot_matrix)
+ composite = N.dot(composite, image_G['mat'])
+
+ # allocate memory from Python as memory leaks when created in C-ext
+ joint_histogram = N.zeros([256, 256], dtype=N.float64);
+
+ if do_lite:
+ R.register_histogram_lite(image_F['data'], image_G['data'], composite, sample_vector, joint_histogram)
+ else:
+ R.register_histogram(image_F['data'], image_G['data'], composite, sample_vector, joint_histogram)
+
+ # smooth the histogram
+ if smooth:
+ p = N.ceil(2*fwhm[0]).astype(int)
+ x = N.array(range(-p, p+1))
+ kernel1 = smooth_kernel(fwhm[0], x)
+ p = N.ceil(2*fwhm[1]).astype(int)
+ x = N.array(range(-p, p+1))
+ kernel2 = smooth_kernel(fwhm[1], x)
+ output=None
+ # 2D filter in 1D separable stages
+ axis = 0
+ result = NDI.correlate1d(joint_histogram, kernel1, axis, output)
+ axis = 1
+ joint_histogram = NDI.correlate1d(result, kernel1, axis, output)
+
+ joint_histogram += epsilon # prevent log(0)
+ # normalize the joint histogram
+ joint_histogram /= joint_histogram.sum()
+ # get the marginals
+ marginal_col = joint_histogram.sum(axis=0)
+ marginal_row = joint_histogram.sum(axis=1)
+
+ if method == 'mi':
+ # mutual information
+ marginal_outer = N.outer(marginal_col, marginal_row)
+ H = joint_histogram * N.log(joint_histogram / marginal_outer)
+ mutual_information = H.sum()
+ cost = -mutual_information
+
+ elif method == 'ecc':
+ # entropy correlation coefficient
+ marginal_outer = N.outer(marginal_col, marginal_row)
+ H = joint_histogram * N.log(joint_histogram / marginal_outer)
+ mutual_information = H.sum()
+ row_entropy = marginal_row * N.log(marginal_row)
+ col_entropy = marginal_col * N.log(marginal_col)
+ ecc = -2.0*mutual_information/(row_entropy.sum() + col_entropy.sum())
+ cost = -ecc
+
+ elif method == 'nmi':
+ # normalized mutual information
+ row_entropy = marginal_row * N.log(marginal_row)
+ col_entropy = marginal_col * N.log(marginal_col)
+ H = joint_histogram * N.log(joint_histogram)
+ nmi = (row_entropy.sum() + col_entropy.sum()) / (H.sum())
+ cost = -nmi
+
+ elif method == 'ncc':
+ # cross correlation from the joint histogram
+ r, c = joint_histogram.shape
+ i = N.array(range(1,c+1))
+ j = N.array(range(1,r+1))
+ m1 = (marginal_row * i).sum()
+ m2 = (marginal_col * j).sum()
+ sig1 = N.sqrt((marginal_row*(N.square(i-m1))).sum())
+ sig2 = N.sqrt((marginal_col*(N.square(j-m2))).sum())
+ [a, b] = N.mgrid[1:c+1, 1:r+1]
+ a = a - m1
+ b = b - m2
+ # element multiplies in the joint histogram and grids
+ H = ((joint_histogram * a) * b).sum()
+ ncc = H / (N.dot(sig1, sig2))
+ cost = -ncc
+
+ if ret_histo:
+ return cost, joint_histogram
+ else:
+ return cost
+
+
+def build_structs(step=2):
+ # build image data structures here
+ P = N.zeros(6, dtype=N.float64);
+ T = N.zeros(6, dtype=N.float64);
+ F = N.zeros(2, dtype=N.int32);
+ S = N.ones(3, dtype=N.int32);
+ sample = N.zeros(2, dtype=N.int32);
+ S[0] = step
+ S[1] = step
+ S[2] = step
+ # histogram smoothing
+ F[0] = 3
+ F[1] = 3
+ # subsample for multiresolution registration
+ sample[0] = 4
+ sample[1] = 2
+ # tolerances for angle (0-2) and translation (3-5)
+ T[0] = 0.02
+ T[1] = 0.02
+ T[2] = 0.02
+ T[3] = 0.001
+ T[4] = 0.001
+ T[5] = 0.001
+ # P[0] = alpha <=> pitch. + alpha is moving back in the sagittal plane
+ # P[1] = beta <=> roll. + beta is moving right in the coronal plane
+ # P[2] = gamma <=> yaw. + gamma is right turn in the transverse plane
+ # P[3] = Tx
+ # P[4] = Ty
+ # P[5] = Tz
+ img_data = {'parms' : P, 'step' : S, 'fwhm' : F, 'tol' : T, 'sample' : sample}
+ return img_data
+
+
+def build_rotate_matrix(img_data_parms):
+ R1 = N.zeros([4,4], dtype=N.float64);
+ R2 = N.zeros([4,4], dtype=N.float64);
+ R3 = N.zeros([4,4], dtype=N.float64);
+ T = N.eye(4, dtype=N.float64);
+
+ alpha = math.radians(img_data_parms[0])
+ beta = math.radians(img_data_parms[1])
+ gamma = math.radians(img_data_parms[2])
+
+ R1[0][0] = 1.0
+ R1[1][1] = math.cos(alpha)
+ R1[1][2] = math.sin(alpha)
+ R1[2][1] = -math.sin(alpha)
+ R1[2][2] = math.cos(alpha)
+ R1[3][3] = 1.0
+
+ R2[0][0] = math.cos(beta)
+ R2[0][2] = math.sin(beta)
+ R2[1][1] = 1.0
+ R2[2][0] = -math.sin(beta)
+ R2[2][2] = math.cos(beta)
+ R2[3][3] = 1.0
+
+ R3[0][0] = math.cos(gamma)
+ R3[0][1] = math.sin(gamma)
+ R3[1][0] = -math.sin(gamma)
+ R3[1][1] = math.cos(gamma)
+ R3[2][2] = 1.0
+ R3[3][3] = 1.0
+
+ T[0][0] = 1.0
+ T[1][1] = 1.0
+ T[2][2] = 1.0
+ T[3][3] = 1.0
+ T[0][3] = img_data_parms[3]
+ T[1][3] = img_data_parms[4]
+ T[2][3] = img_data_parms[5]
+
+ rot_matrix = N.dot(T, R1);
+ rot_matrix = N.dot(rot_matrix, R2);
+ rot_matrix = N.dot(rot_matrix, R3);
+
+ return rot_matrix
+
+
+
+
+
+
Deleted: trunk/scipy/ndimage/registration.py
===================================================================
--- trunk/scipy/ndimage/registration.py 2008-02-13 00:00:26 UTC (rev 3928)
+++ trunk/scipy/ndimage/registration.py 2008-02-13 00:11:57 UTC (rev 3929)
@@ -1,440 +0,0 @@
-import math
-import os
-import numpy as N
-import scipy.ndimage._register as R
-import scipy.special as SP
-import scipy.ndimage as NDI
-import scipy.optimize as OPT
-import time
-
-# anatomical MRI to test with
-inputname = 'ANAT1_V0001.img'
-filename = os.path.join(os.path.split(__file__)[0], inputname)
-
-def remap_image(image, parm_vector):
- M_inverse = get_inverse_mappings(parm_vector)
- # allocate the zero image
- remaped_image = load_blank_image()
- imdata = build_structs(step=1)
- # trilinear interpolation mapping. to be replaced with splines
- R.register_linear_resample(image['data'], remaped_image['data'], M_inverse, imdata['step'])
- return remaped_image
-
-def get_inverse_mappings(parm_vector):
- # get the inverse mapping to rotate the G matrix to F space following registration
- imdata = build_structs(step=1)
- # inverse angles and translations
- imdata['parms'][0] = -parm_vector[0]
- imdata['parms'][1] = -parm_vector[1]
- imdata['parms'][2] = -parm_vector[2]
- imdata['parms'][3] = -parm_vector[3]
- imdata['parms'][4] = -parm_vector[4]
- imdata['parms'][5] = -parm_vector[5]
- M_inverse = build_rotate_matrix(imdata['parms'])
- return M_inverse
-
-def python_coreg(image1, image2, imdata, ftype=1, smimage=0, lite=0, smhist=0,
- method='nmi', opt_method='powell'):
- # image1 is imageF and image2 is imageG in SPM lingo
- # get these from get_test_images for the debug work
- start = time.time()
- # smooth of the images
- if smimage:
- image_F_xyz1 = filter_image_3D(image1['data'], image1['fwhm'], ftype)
- image1['data'] = image_F_xyz1
- image_F_xyz2 = filter_image_3D(image2['data'], image2['fwhm'], ftype)
- image2['data'] = image_F_xyz2
- parm_vector = multires_registration(image1, image2, imdata, lite, smhist, method, opt_method)
- stop = time.time()
- print 'Total Optimizer Time is ', (stop-start)
- return parm_vector
-
-def get_test_images(alpha=0.0, beta=0.0, gamma=0.0):
- image1 = load_image()
- image2 = load_blank_image()
- imdata = build_structs(step=1)
- # allow the G image to be rotated for testing
- imdata['parms'][0] = alpha
- imdata['parms'][1] = beta
- imdata['parms'][2] = gamma
- image1['fwhm'] = build_fwhm(image1['mat'], imdata['step'])
- image2['fwhm'] = build_fwhm(image2['mat'], imdata['step'])
- M = build_rotate_matrix(imdata['parms'])
- R.register_linear_resample(image1['data'], image2['data'], M, imdata['step'])
- return image1, image2, imdata
-
-def multires_registration(image1, image2, imdata, lite, smhist, method, opt_method):
- ret_histo=0
- # zero out the start parameter; but this may be set to large values
- # if the head is out of range and well off the optimal alignment skirt
- imdata['parms'][0:5] = 0.0
- # make the step a scalar to can put in a multi-res loop
- loop = range(imdata['sample'].size)
- x = imdata['parms']
- for i in loop:
- step = imdata['sample'][i]
- imdata['step'][:] = step
- optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
- p_args = (optfunc_args,)
- if opt_method=='powell':
- print 'POWELL multi-res registration step size ', step
- print 'vector ', x
- x = OPT.fmin_powell(optimize_function, x, args=p_args, callback=callback_powell)
- elif opt_method=='cg':
- print 'CG multi-res registration step size ', step
- print 'vector ', x
- x = OPT.fmin_cg(optimize_function, x, args=p_args, callback=callback_cg)
- elif opt_method=='hybrid':
- if i==0:
- print 'Hybrid POWELL multi-res registration step size ', step
- print 'vector ', x
- lite = 0
- optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
- p_args = (optfunc_args,)
- x = OPT.fmin_powell(optimize_function, x, args=p_args, callback=callback_powell)
- elif i==1:
- print 'Hybrid CG multi-res registration step size ', step
- print 'vector ', x
- lite = 1
- optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
- p_args = (optfunc_args,)
- x = OPT.fmin_cg(optimize_function, x, args=p_args, callback=callback_cg)
-
- return x
-
-
-def test_image_filter(image, imdata, ftype=2):
- # test the 3D image filter on an image. ftype 1 is SPM, ftype 2 is simple Gaussian
- image['fwhm'] = build_fwhm(image['mat'], imdata['step'])
- filt_image = filter_image_3D(image['data'], image['fwhm'], ftype)
- return filt_image
-
-def callback_powell(x):
- print 'Parameter Vector from Powell: - '
- print x
- return
-
-def callback_cg(x):
- print 'Parameter Vector from Conjugate Gradient: - '
- print x
- return
-
-def test_alignment(image1, image2, imdata, method='ncc', lite=0, smhist=0,
- alpha=0.0, beta=0.0, gamma=0.0, ret_histo=0):
-
- # to test the cost function and view the joint histogram
- # for 2 images. used for debug
- imdata['parms'][0] = alpha
- imdata['parms'][1] = beta
- imdata['parms'][2] = gamma
- M = build_rotate_matrix(imdata['parms'])
- optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
-
- if ret_histo:
- cost, joint_histogram = optimize_function(imdata['parms'], optfunc_args)
- return cost, joint_histogram
- else:
- cost = optimize_function(imdata['parms'], optfunc_args)
- return cost
-
-
-def smooth_kernel(fwhm, x, ktype=1):
- eps = 0.00001
- s = N.square((fwhm/math.sqrt(8.0*math.log(2.0)))) + eps
- if ktype==1:
- # from SPM: Gauss kernel convolved with 1st degree B spline
- w1 = 0.5 * math.sqrt(2.0/s)
- w2 = -0.5 / s
- w3 = math.sqrt((s*math.pi) /2.0)
- kernel = 0.5*(SP.erf(w1*(x+1))*(x+1) + SP.erf(w1*(x-1))*(x-1) - 2.0*SP.erf(w1*(x))*(x) +
- w3*(N.exp(w2*N.square(x+1))) + N.exp(w2*(N.square(x-1))) - 2.0*N.exp(w2*N.square(x)))
- kernel[kernel<0] = 0
- kernel = kernel / kernel.sum()
- else:
- # Gauss kernel
- kernel = (1.0/math.sqrt(2.0*math.pi*s)) * N.exp(-N.square(x)/(2.0*s))
- kernel = kernel / kernel.sum()
-
- return kernel
-
-def filter_image_3D(imageRaw, fwhm, ftype=2):
- p = N.ceil(2*fwhm[0]).astype(int)
- x = N.array(range(-p, p+1))
- kernel_x = smooth_kernel(fwhm[0], x, ktype=ftype)
- p = N.ceil(2*fwhm[1]).astype(int)
- x = N.array(range(-p, p+1))
- kernel_y = smooth_kernel(fwhm[1], x, ktype=ftype)
- p = N.ceil(2*fwhm[2]).astype(int)
- x = N.array(range(-p, p+1))
- kernel_z = smooth_kernel(fwhm[2], x, ktype=ftype)
- output=None
- # 3D filter in 3 1D separable stages
- axis = 0
- image_F_x = NDI.correlate1d(imageRaw, kernel_x, axis, output)
- axis = 1
- image_F_xy = NDI.correlate1d(image_F_x, kernel_y, axis, output)
- axis = 2
- image_F_xyz = NDI.correlate1d(image_F_xy, kernel_z, axis, output)
- return image_F_xyz
-
-
-def resample_image(smimage=0, ftype=2, alpha=0.0, beta=0.0, gamma=0.0,
- Tx=0.0, Ty=0.0, Tz=0.0, stepsize=1):
-
- # takes an image and 3D rotate using trilinear interpolation
- image1 = load_image()
- image2 = load_blank_image()
- imdata = build_structs(step=stepsize)
- imdata['parms'][0] = alpha
- imdata['parms'][1] = beta
- imdata['parms'][2] = gamma
- imdata['parms'][3] = Tx
- imdata['parms'][4] = Ty
- imdata['parms'][5] = Tz
- image1['fwhm'] = build_fwhm(image1['mat'], imdata['step'])
- image2['fwhm'] = build_fwhm(image2['mat'], imdata['step'])
- M = build_rotate_matrix(imdata['parms'])
- if smimage:
- image_F_xyz1 = filter_image_3D(image1['data'], image1['fwhm'], ftype)
- image1['data'] = image_F_xyz1
- image_F_xyz2 = filter_image_3D(image2['data'], image2['fwhm'], ftype)
- image2['data'] = image_F_xyz2
-
- R.register_linear_resample(image1['data'], image2['data'], M, imdata['step'])
-
- return image2
-
-
-def build_fwhm(M, S):
- view_3x3 = N.square(M[0:3, 0:3])
- vxg = N.sqrt(view_3x3.sum(axis=0))
- # assumes that sampling is the same for xyz
- size = N.array([1,1,1])*S[0]
- x = N.square(size) - N.square(vxg)
- # clip
- x[x<0] = 0
- fwhm = N.sqrt(x) / vxg
- # pathology when stepsize = 1 for MAT equal to the identity matrix
- fwhm[fwhm==0] = 1
- # return the 3D Gaussian kernel width (xyz)
- return fwhm
-
-def load_image(imagename=filename, rows=256, cols=256, layers=90):
- ImageVolume = N.fromfile(imagename, dtype=N.uint16).reshape(layers, rows, cols);
- # clip to 8 bits. this will be rescale to 8 bits for fMRI
- ImageVolume[ImageVolume>255] = 255
- # voxel to pixel is identity for this simulation using anatomical MRI volume
- # 4x4 matrix
- M = N.eye(4, dtype=N.float64);
- # dimensions
- D = N.zeros(3, dtype=N.int32);
- # Gaussian kernel - fill in with build_fwhm()
- F = N.zeros(3, dtype=N.float64);
- D[0] = rows
- D[1] = cols
- D[2] = layers
- # make sure the data type is uchar
- ImageVolume = ImageVolume.astype(N.uint8)
- image = {'data' : ImageVolume, 'mat' : M, 'dim' : D, 'fwhm' : F}
- return image
-
-def load_blank_image(rows=256, cols=256, layers=90):
- ImageVolume = N.zeros(layers*rows*cols, dtype=N.uint16).reshape(layers, rows, cols);
- # voxel to pixel is identity for this simulation using anatomical MRI volume
- # 4x4 matrix
- M = N.eye(4, dtype=N.float64);
- # dimensions
- D = N.zeros(3, dtype=N.int32);
- # Gaussian kernel - fill in with build_fwhm()
- F = N.zeros(3, dtype=N.float64);
- D[0] = rows
- D[1] = cols
- D[2] = layers
- # make sure the data type is uchar
- ImageVolume = ImageVolume.astype(N.uint8)
- image = {'data' : ImageVolume, 'mat' : M, 'dim' : D, 'fwhm' : F}
- return image
-
-def optimize_function(x, optfunc_args):
- image_F = optfunc_args[0]
- image_G = optfunc_args[1]
- sample_vector = optfunc_args[2]
- fwhm = optfunc_args[3]
- do_lite = optfunc_args[4]
- smooth = optfunc_args[5]
- method = optfunc_args[6]
- ret_histo = optfunc_args[7]
-
- rot_matrix = build_rotate_matrix(x)
- cost = 0.0
- epsilon = 2.2e-16
- # image_G is base image
- # image_F is the to-be-rotated image
- # rot_matrix is the 4x4 constructed (current angles and translates) transform matrix
- # sample_vector is the subsample vector for x-y-z
-
- F_inv = N.linalg.inv(image_F['mat'])
- composite = N.dot(F_inv, rot_matrix)
- composite = N.dot(composite, image_G['mat'])
-
- # allocate memory from Python as memory leaks when created in C-ext
- joint_histogram = N.zeros([256, 256], dtype=N.float64);
-
- if do_lite:
- R.register_histogram_lite(image_F['data'], image_G['data'], composite, sample_vector, joint_histogram)
- else:
- R.register_histogram(image_F['data'], image_G['data'], composite, sample_vector, joint_histogram)
-
- # smooth the histogram
- if smooth:
- p = N.ceil(2*fwhm[0]).astype(int)
- x = N.array(range(-p, p+1))
- kernel1 = smooth_kernel(fwhm[0], x)
- p = N.ceil(2*fwhm[1]).astype(int)
- x = N.array(range(-p, p+1))
- kernel2 = smooth_kernel(fwhm[1], x)
- output=None
- # 2D filter in 1D separable stages
- axis = 0
- result = NDI.correlate1d(joint_histogram, kernel1, axis, output)
- axis = 1
- joint_histogram = NDI.correlate1d(result, kernel1, axis, output)
-
- joint_histogram += epsilon # prevent log(0)
- # normalize the joint histogram
- joint_histogram /= joint_histogram.sum()
- # get the marginals
- marginal_col = joint_histogram.sum(axis=0)
- marginal_row = joint_histogram.sum(axis=1)
-
- if method == 'mi':
- # mutual information
- marginal_outer = N.outer(marginal_col, marginal_row)
- H = joint_histogram * N.log(joint_histogram / marginal_outer)
- mutual_information = H.sum()
- cost = -mutual_information
-
- elif method == 'ecc':
- # entropy correlation coefficient
- marginal_outer = N.outer(marginal_col, marginal_row)
- H = joint_histogram * N.log(joint_histogram / marginal_outer)
- mutual_information = H.sum()
- row_entropy = marginal_row * N.log(marginal_row)
- col_entropy = marginal_col * N.log(marginal_col)
- ecc = -2.0*mutual_information/(row_entropy.sum() + col_entropy.sum())
- cost = -ecc
-
- elif method == 'nmi':
- # normalized mutual information
- row_entropy = marginal_row * N.log(marginal_row)
- col_entropy = marginal_col * N.log(marginal_col)
- H = joint_histogram * N.log(joint_histogram)
- nmi = (row_entropy.sum() + col_entropy.sum()) / (H.sum())
- cost = -nmi
-
- elif method == 'ncc':
- # cross correlation from the joint histogram
- r, c = joint_histogram.shape
- i = N.array(range(1,c+1))
- j = N.array(range(1,r+1))
- m1 = (marginal_row * i).sum()
- m2 = (marginal_col * j).sum()
- sig1 = N.sqrt((marginal_row*(N.square(i-m1))).sum())
- sig2 = N.sqrt((marginal_col*(N.square(j-m2))).sum())
- [a, b] = N.mgrid[1:c+1, 1:r+1]
- a = a - m1
- b = b - m2
- # element multiplies in the joint histogram and grids
- H = ((joint_histogram * a) * b).sum()
- ncc = H / (N.dot(sig1, sig2))
- cost = -ncc
-
- if ret_histo:
- return cost, joint_histogram
- else:
- return cost
-
-
-def build_structs(step=2):
- # build image data structures here
- P = N.zeros(6, dtype=N.float64);
- T = N.zeros(6, dtype=N.float64);
- F = N.zeros(2, dtype=N.int32);
- S = N.ones(3, dtype=N.int32);
- sample = N.zeros(2, dtype=N.int32);
- S[0] = step
- S[1] = step
- S[2] = step
- # histogram smoothing
- F[0] = 3
- F[1] = 3
- # subsample for multiresolution registration
- sample[0] = 4
- sample[1] = 2
- # tolerances for angle (0-2) and translation (3-5)
- T[0] = 0.02
- T[1] = 0.02
- T[2] = 0.02
- T[3] = 0.001
- T[4] = 0.001
- T[5] = 0.001
- # P[0] = alpha <=> pitch. + alpha is moving back in the sagittal plane
- # P[1] = beta <=> roll. + beta is moving right in the coronal plane
- # P[2] = gamma <=> yaw. + gamma is right turn in the transverse plane
- # P[3] = Tx
- # P[4] = Ty
- # P[5] = Tz
- img_data = {'parms' : P, 'step' : S, 'fwhm' : F, 'tol' : T, 'sample' : sample}
- return img_data
-
-
-def build_rotate_matrix(img_data_parms):
- R1 = N.zeros([4,4], dtype=N.float64);
- R2 = N.zeros([4,4], dtype=N.float64);
- R3 = N.zeros([4,4], dtype=N.float64);
- T = N.eye(4, dtype=N.float64);
-
- alpha = math.radians(img_data_parms[0])
- beta = math.radians(img_data_parms[1])
- gamma = math.radians(img_data_parms[2])
-
- R1[0][0] = 1.0
- R1[1][1] = math.cos(alpha)
- R1[1][2] = math.sin(alpha)
- R1[2][1] = -math.sin(alpha)
- R1[2][2] = math.cos(alpha)
- R1[3][3] = 1.0
-
- R2[0][0] = math.cos(beta)
- R2[0][2] = math.sin(beta)
- R2[1][1] = 1.0
- R2[2][0] = -math.sin(beta)
- R2[2][2] = math.cos(beta)
- R2[3][3] = 1.0
-
- R3[0][0] = math.cos(gamma)
- R3[0][1] = math.sin(gamma)
- R3[1][0] = -math.sin(gamma)
- R3[1][1] = math.cos(gamma)
- R3[2][2] = 1.0
- R3[3][3] = 1.0
-
- T[0][0] = 1.0
- T[1][1] = 1.0
- T[2][2] = 1.0
- T[3][3] = 1.0
- T[0][3] = img_data_parms[3]
- T[1][3] = img_data_parms[4]
- T[2][3] = img_data_parms[5]
-
- rot_matrix = N.dot(T, R1);
- rot_matrix = N.dot(rot_matrix, R2);
- rot_matrix = N.dot(rot_matrix, R3);
-
- return rot_matrix
-
-
-
-
-
-
More information about the Scipy-svn
mailing list