[Scipy-svn] r3969 - trunk/scipy/ndimage

scipy-svn at scipy.org scipy-svn at scipy.org
Tue Mar 4 21:14:01 EST 2008


Author: tom.waite
Date: 2008-03-04 20:13:58 -0600 (Tue, 04 Mar 2008)
New Revision: 3969

Modified:
   trunk/scipy/ndimage/_registration.py
Log:
Added fMRI series co-registration and resampling.

Modified: trunk/scipy/ndimage/_registration.py
===================================================================
--- trunk/scipy/ndimage/_registration.py	2008-03-04 03:58:32 UTC (rev 3968)
+++ trunk/scipy/ndimage/_registration.py	2008-03-05 02:13:58 UTC (rev 3969)
@@ -75,7 +75,7 @@
 	remaped_image = N.zeros(layers*rows*cols, dtype=N.uint8).reshape(layers, rows, cols)
 	remaped_image = {'data' : remaped_image, 'mat' : image['mat'], 
 			 'dim' : image['dim'], 'fwhm' : image['fwhm']}
-	imdata = build_structs(step=1)
+	imdata = build_structs()
 
 	if resample == 'linear':
 	    # trilinear interpolation mapping.
@@ -88,7 +88,7 @@
 
 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)
+	imdata = build_structs()
 	# inverse angles and translations
 	imdata['parms'][0] = -parm_vector[0]
 	imdata['parms'][1] = -parm_vector[1]
@@ -106,8 +106,6 @@
     	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)
@@ -210,13 +208,13 @@
 
 
 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): 
+		   Tx=0.0, Ty=0.0, Tz=0.0): 
 	            
 	# takes an image and 3D rotate using trilinear interpolation
 	anat_desc = load_anatMRI_desc()
 	image1 = load_volume(anat_desc, imagename='ANAT1_V0001.img')
 	image2 = load_volume(anat_desc, imagename=None)
-	imdata = build_structs(step=stepsize)
+	imdata = build_structs()
 	imdata['parms'][0] = alpha
 	imdata['parms'][1] = beta
 	imdata['parms'][2] = gamma
@@ -270,13 +268,9 @@
 	# 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'])
 	F_inv = N.linalg.inv(image_F['mat'])
 	composite = N.dot(F_inv, image_G['mat'])
 	composite = N.dot(composite, rot_matrix)
-	#print ' composite ', composite 
 
 	# allocate memory from Python as memory leaks when created in C-ext
 	joint_histogram = N.zeros([256, 256], dtype=N.float64);
@@ -436,8 +430,10 @@
 
 
 def load_volume(imagedesc, imagename=None, threshold=0.999, debug=0):
-	# imagename of none means to create a blank image
+	# load MRI or fMRI volume and return an autoscaled 8 bit image.
+	# autoscale is using integrated histogram to deal with outlier high amplitude voxels
 	if imagename == None:
+	    # imagename of none means to create a blank image
     	    ImageVolume = N.zeros(imagedesc['layers']*imagedesc['rows']*imagedesc['cols'],
 			   dtype=N.uint16).reshape(imagedesc['layers'], imagedesc['rows'], imagedesc['cols'])
 	else:
@@ -514,33 +510,40 @@
 	files_fMRI = glob.glob(path)
 	return files_fMRI
 
-def build_aligned_fMRI_mean_volume(path):
-	desc = load_fMRI_desc()
-	ave_fMRI_volume = N.zeros(desc['layers']*desc['rows']*desc['cols'],
-			  dtype=N.float64).reshape(desc['layers'], desc['rows'], desc['cols'])
-	data = read_fMRI_directory(path)
-	count = 0
-	for i in data:
-	    print 'add volume ', i
-	    # this uses integrated histogram normalization
-	    image = load_volume(desc, i)
-	    # co-reg successive pairs, then remap and ave
-	    ave_fMRI_volume = ave_fMRI_volume + image['data'].astype(N.float64)
-	    count = count + 1
 
-	ave_fMRI_volume = ave_fMRI_volume / float(count)
+def check_alignment(image1, image2, imdata, method='ncc', lite=0, smhist=0, 
+		    alpha=0.0, beta=0.0, gamma=0.0, Tx=0, Ty=0, Tz=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
+	M = build_rotate_matrix(imdata['parms'])
+	optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
 
-	return ave_fMRI_volume 
+	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
 
+
 #
-#  ---- testing/debug routines  ---- 
+#  ---- demo/debug routines  ---- 
 #
 
 def build_scale_image(image, scale):
-	(layers, rows, cols) = image['data'].shape
 	#
 	# rescale the 'mat' (voxel to physical mapping matrix) 
 	#
+	(layers, rows, cols) = image['data'].shape
 	M = image['mat'] * scale
 	# dimensions 
 	D = N.zeros(3, dtype=N.int32);
@@ -557,7 +560,7 @@
 	return scaled_image
 
 
-def get_test_MRI_volumes(scale=2, alpha=3.0, beta=4.0, gamma=5.0):
+def demo_MRI_volume_align(scale=2, alpha=3.0, beta=4.0, gamma=5.0, Tx = 0.0, Ty = 0.0, Tz = 0.0):
 	#
 	# this is for coreg MRI / fMRI scale test. The volume is anatomical MRI.
 	# the image is rotated in 3D. after rotation the image is scaled.  
@@ -566,64 +569,142 @@
 	anat_desc = load_anatMRI_desc()
 	image1 = load_volume(anat_desc, imagename='ANAT1_V0001.img')
 	image2 = load_volume(anat_desc, imagename=None)
-	imdata = build_structs(step=1)
+	imdata = build_structs()
 	image1['fwhm'] = build_fwhm(image1['mat'], imdata['step'])
 	image2['fwhm'] = build_fwhm(image2['mat'], imdata['step'])
 	imdata['parms'][0] = alpha
 	imdata['parms'][1] = beta
 	imdata['parms'][2] = gamma
+	imdata['parms'][3] = Tx
+	imdata['parms'][4] = Ty
+	imdata['parms'][5] = Tz
 	M = build_rotate_matrix(imdata['parms'])
+	# rotate volume. linear interpolation means the volume is low pass filtered
 	R.register_linear_resample(image1['data'], image2['data'], M, imdata['step'])
+	# subsample volume
 	image3 = build_scale_image(image2, scale)
 	return image1, image3, imdata
 
-def get_test_fMRI_rotated(fMRIVol, alpha=3.0, beta=4.0, gamma=5.0):
+def demo_rotate_fMRI_volume(fMRIVol, x): 
 	#
-	# return rotated fMRIVol
+	# return rotated fMRIVol. the fMRIVol is already loaded, and gets rotated
 	#
 
 	desc = load_fMRI_desc()
 	image = load_volume(desc, imagename=None)
-	imdata = build_structs(step=1)
+	imdata = build_structs()
 	image['fwhm'] = build_fwhm(image['mat'], imdata['step'])
-	imdata['parms'][0] = alpha
-	imdata['parms'][1] = beta
-	imdata['parms'][2] = gamma
+	imdata['parms'][0] = x[0]  # alpha
+	imdata['parms'][1] = x[1]  # beta
+	imdata['parms'][2] = x[2]  # gamma
+	imdata['parms'][3] = x[3]  # Tx
+	imdata['parms'][4] = x[4]  # Ty
+	imdata['parms'][5] = x[5]  # Tz
 	M = build_rotate_matrix(imdata['parms'])
+	# rotate volume. cubic spline interpolation means the volume is NOT low pass filtered
 	R.register_cubic_resample(fMRIVol['data'], image['data'], M, imdata['step'])
 	return image
 
+def demo_MRI_coregistration(optimizer_method='powell', histo_method=1, smooth_histo=0, smooth_image=0, ftype=1):
+	# demo of alignment of fMRI series with anatomical MRI
+	# in this demo, each fMRI volume is first perturbed (rotated, translated) 
+	# by a random value. The initial registration is measured, then the optimal
+	# alignment is computed and the registration measure made following the volume remap.
+	# The fMRI registration is done with the first fMRI volume using normalized cross-correlation.
+	# Each fMRI volume is rotated to the fMRI-0 volume and the series is ensemble averaged.
+	# The ensemble averaged is then registered with the anatomical MRI volume using normalized mutual information.
+	# The fMRI series is then rotated with this parameter. The alignments are done with 3D cubic splines.
 
-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
+	# read the anatomical MRI volume
+	anat_desc = load_anatMRI_desc()
+	imageF_anat = load_volume(anat_desc, imagename='ANAT1_V0001.img')
+	# the sampling structure
+	imdata = build_structs()
+	# the volume filter
+	imageF_anat['fwhm'] = build_fwhm(imageF_anat['mat'], imdata['step'])
 
+	# read in the file list of the fMRI data
+	metric_test = N.dtype([('cost', 'f'),
+                    	       ('align_cost', 'f'),
+                    	       ('rotate', 'f', 6),
+                    	       ('align_rotate', 'f', 6)])
 
-def test_alignment(image1, image2, imdata, method='ncc', lite=0, smhist=0, 
-		        alpha=0.0, beta=0.0, gamma=0.0, Tx=0, Ty=0, Tz=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
-	M = build_rotate_matrix(imdata['parms'])
-	optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
+	fMRIdata = read_fMRI_directory('fMRIData\*.img')
+	fmri_desc = load_fMRI_desc()
+	fmri_series = {}
+	ave_fMRI_volume = N.zeros(fmri_desc['layers']*fmri_desc['rows']*fmri_desc['cols'],
+			  dtype=N.float64).reshape(fmri_desc['layers'], fmri_desc['rows'], fmri_desc['cols'])
+	count = 0
+	number_volumes = len(fMRIdata)
+	measures = N.zeros(number_volumes, dtype=metric_test)
+	# load and perturb (rotation, translation) the fMRI volumes
+	for i in fMRIdata:
+	    image = load_volume(fmri_desc, i)
+	    # random perturbation of angle, translation for each volume beyond the first
+	    if count == 0:
+		image['fwhm'] = build_fwhm(image['mat'], imdata['step'])
+	        fmri_series[count] = image
+	        count = count + 1
+	    else:
+	        x = N.random.random(6) - 0.5
+	        x = 10.0 * x
+	        fmri_series[count] = demo_rotate_fMRI_volume(image, x)
+		measures[count]['rotate'][0:6] = x[0:6]
+	        count = count + 1
 
-	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
 
+	# load and register the fMRI volumes with volume_0 using normalized cross correlation metric
+	imageF = fmri_series[0]
+	if smooth_image:
+	    image_F_xyz = filter_image_3D(imageF['data'], imageF['fwhm'], ftype)
+	    imageF['data'] = image_F_xyz
+	for i in range(1, number_volumes):
+	    imageG = fmri_series[i]
+	    # the measure prior to alignment 
+	    measures[i]['cost'] = check_alignment(imageF, imageG, imdata, method='ncc',
+			                          lite=histo_method, smhist=smooth_histo)
+	    x = python_coreg(imageF, imageG, imdata, lite=histo_method, method='ncc',
+			     opt_method=optimizer_method, smhist=smooth_histo, smimage=smooth_image)
+	    measures[i]['align_rotate'][0:6] = x[0:6]
+	    measures[i]['align_cost'] = check_alignment(imageF, imageG, imdata, method='ncc', 
+		                             lite=histo_method, smhist=smooth_histo,
+					     alpha=x[0], beta=x[1], gamma=x[2], Tx=x[3], Ty=x[4], Tz=x[5])
 
+
+	# align the volumes and average them for co-registration with the anatomical MRI 
+	ave_fMRI_volume = fmri_series[0]['data'].astype(N.float64)
+	for i in range(1, number_volumes):
+	    image = fmri_series[i]
+	    x[0:6] = measures[i]['align_rotate'][0:6]
+	    # overwrite the fMRI volume with the aligned volume
+	    fmri_series[i] = remap_image(image, x, resample='cubic')
+	    ave_fMRI_volume = ave_fMRI_volume + fmri_series[i]['data'].astype(N.float64)
+
+	ave_fMRI_volume = (ave_fMRI_volume / float(number_volumes)).astype(N.uint8)
+	ave_fMRI_volume = {'data' : ave_fMRI_volume, 'mat' : imageF['mat'], 
+			   'dim' : imageF['dim'], 'fwhm' : imageF['fwhm']}
+	# register (using normalized mutual information) with the anatomical MRI
+	if smooth_image:
+	    image_F_anat_xyz = filter_image_3D(imageF_anat['data'], imageF_anat['fwhm'], ftype)
+	    imageF_anat['data'] = image_F_anat_xyz
+	x = python_coreg(imageF_anat, ave_fMRI_volume, imdata, lite=histo_method,
+			 method='nmi', opt_method=optimizer_method, smhist=smooth_histo, smimage=smooth_image)
+	print 'functional-anatomical align parameters '
+	print x
+	for i in range(number_volumes):
+	    image = fmri_series[i]
+	    # overwrite the fMRI volume with the anatomical-aligned volume
+	    fmri_series[i] = remap_image(image, x, resample='cubic')
+
+	return measures, imageF_anat, fmri_series
+
+
+def demo_fMRI_resample(imageF_anat, fmri_series):
+	resampled_fmri_series = {}
+	number_volumes = len(fmri_series)
+	for i in range(number_volumes):
+	    resampled_fmri_series[i] = resize_image(fmri_series[i], imageF_anat['mat'])
+
+	return resampled_fmri_series
+
+




More information about the Scipy-svn mailing list