[Scipy-svn] r3896 - trunk/scipy/ndimage
scipy-svn@scip...
scipy-svn@scip...
Tue Feb 5 19:18:30 CST 2008
Author: tom.waite
Date: 2008-02-05 19:18:26 -0600 (Tue, 05 Feb 2008)
New Revision: 3896
Modified:
trunk/scipy/ndimage/registration.py
Log:
New functionality for NIPY registration.
Modified: trunk/scipy/ndimage/registration.py
===================================================================
--- trunk/scipy/ndimage/registration.py 2008-02-06 01:17:39 UTC (rev 3895)
+++ trunk/scipy/ndimage/registration.py 2008-02-06 01:18:26 UTC (rev 3896)
@@ -8,36 +8,50 @@
import time
# anatomical MRI to test with
-# test registration on same image (optimal vector is (0,0,0,0,0,0)
inputname = 'ANAT1_V0001.img'
filename = os.path.join(os.path.split(__file__)[0], inputname)
-def python_coreg(ftype=2, smimage=0, lite=1, smhist=0, method='mi', opt_method='powell'):
- # get_images is testing with 2 copies of anatomical MRI
- image1, image2, imdata = get_images(ftype, smimage)
+def get_mappings(parm_vector):
+ # get the inverse mapping to rotate the G matrix to F space following registration
+ M_foreward = build_rotate_matrix(parm_vector)
+ M_inverse = N.linalg.inv(M_foreward)
+ return M_foreward, 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_images(ftype, smimage):
+def get_test_images(alpha=0.0, beta=0.0, gamma=0.0):
image1 = load_image()
- image2 = load_image()
- imdata = build_structs()
+ 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'])
- 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
-
+ 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']
@@ -54,35 +68,31 @@
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 test_optimizers(step=2, smooth=0, shist=0):
- opt_stats = {}
- print 'powell with stochastic resampling'
- x_0, p_time_0 = optimizer_powell(lite=0, smimage=smooth, smhist=shist, stepsize=step)
- opt_stats[0] = {'parms' : x_0, 'time' : p_time_0,
- 'label' : 'powell with stochastic resampling'}
- print 'powell without stochastic resampling'
- x_1, p_time_1 = optimizer_powell(lite=1, smimage=smooth, smhist=shist, stepsize=step)
- opt_stats[1] = {'parms' : x_1, 'time' : p_time_1,
- 'label' : 'powell without stochastic resampling'}
- print 'conjugate gradient with stochastic resampling'
- xcg_0, cg_time_0 = optimizer_cg(lite=0, smimage=smooth, smhist=shist, stepsize=step)
- opt_stats[2] = {'parms' : xcg_0, 'time' : cg_time_0,
- 'label' : 'conjugate gradient with stochastic resampling'}
- print 'conjugate gradient without stochastic resampling'
- xcg_1, cg_time_1 = optimizer_cg(lite=1, smimage=smooth, smhist=shist, stepsize=step)
- opt_stats[3] = {'parms' : xcg_1, 'time' : cg_time_1,
- 'label' : 'conjugate gradient without stochastic resampling'}
- return opt_stats
-
def callback_powell(x):
print 'Parameter Vector from Powell: - '
print x
@@ -93,74 +103,15 @@
print x
return
-def optimizer_powell(lite=0, smhist=0, smimage=1, method='mi', ftype=2, stepsize=2):
- # test the Powell registration on the anatomical MRI volume
- image1 = load_image()
- image2 = load_image()
- imdata = build_structs(step=stepsize)
- 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
-
- ret_histo=0
- optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
- p_args = (optfunc_args,)
- start = time.time()
- x = OPT.fmin_powell(optimize_function, imdata['parms'], args=p_args, callback=callback_powell)
- stop = time.time()
- return x, (stop-start)
-
-
-def optimizer_cg(lite=0, smhist=0, smimage=1, method='mi', ftype=2, stepsize=2):
- # test the CG registration on the anatomical MRI volume
- image1 = load_image()
- image2 = load_image()
- imdata = build_structs(step=stepsize)
- 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
-
- ret_histo=0
- optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
- p_args = (optfunc_args,)
- start = time.time()
- x = OPT.fmin_cg(optimize_function, imdata['parms'], args=p_args, callback=callback_cg)
- stop = time.time()
- return x, (stop-start)
-
-
-def reg_single_pass(lite=0, smhist=0, smimage=0, method='mi', ftype=2, alpha=0.0,
- beta=0.0, gamma=0.0, Tx=0.0, Ty=0.0, Tz=0.0, ret_histo=0, stepsize=2):
- image1 = load_image()
- image2 = load_image()
- imdata = build_structs(step=stepsize)
+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
- 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'])
- print 'image1[fwhm] ', image1['fwhm']
- print 'image2[fwhm] ', image2['fwhm']
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
-
optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
if ret_histo:
@@ -210,6 +161,34 @@
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))
@@ -243,6 +222,22 @@
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]
@@ -257,8 +252,8 @@
rot_matrix = build_rotate_matrix(x)
cost = 0.0
epsilon = 2.2e-16
- # image_F is base image
- # image_G is the rotated image
+ # 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
@@ -274,6 +269,7 @@
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))
@@ -366,9 +362,9 @@
T[3] = 0.001
T[4] = 0.001
T[5] = 0.001
- # P[0] = alpha <=> pitch
- # P[1] = beta <=> roll
- # P[2] = gamma <=> yaw
+ # 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
More information about the Scipy-svn
mailing list