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

scipy-svn@scip... scipy-svn@scip...
Fri Feb 29 18:44:37 CST 2008


Author: tom.waite
Date: 2008-02-29 18:44:34 -0600 (Fri, 29 Feb 2008)
New Revision: 3966

Modified:
   trunk/scipy/ndimage/_registration.py
Log:
Bug fix and enhancements

Modified: trunk/scipy/ndimage/_registration.py
===================================================================
--- trunk/scipy/ndimage/_registration.py	2008-03-01 00:44:16 UTC (rev 3965)
+++ trunk/scipy/ndimage/_registration.py	2008-03-01 00:44:34 UTC (rev 3966)
@@ -29,10 +29,52 @@
 inputname = 'ANAT1_V0001.img'
 filename = os.path.join(os.path.split(__file__)[0], inputname)
 
+#
+#  ---- co-registration and IO  ---- 
+#
+
+def resize_image(imageG, imageF_mat):
+	#
+	# Fractional resample imageG to imageF size. 
+	#
+	Z = N.zeros(3, dtype=N.float64);
+	# get the zoom
+	Z[0] = imageG['mat'][0][0] / imageF_mat[0][0]
+	Z[1] = imageG['mat'][1][1] / imageF_mat[1][1]
+	Z[2] = imageG['mat'][2][2] / imageF_mat[2][2]
+
+	# new volume dimensions (rounded)
+	D = N.zeros(3, dtype=N.int32);
+	D[0] = int(float(imageG['dim'][0])*Z[0]+0.5)
+	D[1] = int(float(imageG['dim'][1])*Z[1]+0.5)
+	D[2] = int(float(imageG['dim'][2])*Z[2]+0.5)
+
+	M = N.eye(4, dtype=N.float64);
+	# for the test data, set the xyz voxel sizes for fMRI volume
+	M[0][0] = imageG['mat'][0][0]/Z[0]
+	M[1][1] = imageG['mat'][1][1]/Z[1]
+	M[2][2] = imageG['mat'][2][2]/Z[2]
+
+    	image = N.zeros(D[2]*D[1]*D[0], dtype=N.uint8).reshape(D[2], D[0], D[1])
+	mode  = 2
+	scale = 0
+	R.register_volume_resample(imageG['data'], image, Z, scale, mode)
+	F = N.zeros(3, dtype=N.float64);
+	zoom_image = {'data' : image, 'mat' : M, 'dim' : D, 'fwhm' : F}
+
+	return zoom_image
+
 def remap_image(image, parm_vector, resample='linear'):
+	#
+	# remap imageG to coordinates of imageF (creates imageG')
+	# use the 6 dim parm_vector (3 angles, 3 translations) to remap
+	#
 	M_inverse = get_inverse_mappings(parm_vector)
+	(layers, rows, cols) = image['data'].shape
 	# allocate the zero image
-	remaped_image = load_blank_image()
+	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)
 
 	if resample == 'linear':
@@ -117,12 +159,6 @@
 	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
@@ -133,25 +169,6 @@
 	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
@@ -196,8 +213,9 @@
 		   Tx=0.0, Ty=0.0, Tz=0.0, stepsize=1): 
 	            
 	# takes an image and 3D rotate using trilinear interpolation
-	image1 = load_anatMRI_image()
-	image2 = load_blank_image()
+	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['parms'][0] = alpha
 	imdata['parms'][1] = beta
@@ -252,9 +270,13 @@
 	# 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, rot_matrix)
-	composite = N.dot(composite, image_G['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);
@@ -334,7 +356,7 @@
 	    return cost
 
 
-def build_structs(step=2):
+def build_structs(step=1):
 	# build image data structures here
 	P = N.zeros(6, dtype=N.float64);
 	T = N.zeros(6, dtype=N.float64);
@@ -413,59 +435,36 @@
 	return rot_matrix
 
 
-def load_fMRI_image(imagename, rows=64, cols=64, layers=28, threshold=0.999, debug=0):
-	# un-scaled images
-    	ImageVolume = N.fromfile(imagename, dtype=N.uint16).reshape(layers, rows, cols);
+def load_volume(imagedesc, imagename=None, threshold=0.999, debug=0):
+	# imagename of none means to create a blank image
+	if imagename == None:
+    	    ImageVolume = N.zeros(imagedesc['layers']*imagedesc['rows']*imagedesc['cols'],
+			   dtype=N.uint16).reshape(imagedesc['layers'], imagedesc['rows'], imagedesc['cols'])
+	else:
+    	    ImageVolume = N.fromfile(imagename,
+			   dtype=N.uint16).reshape(imagedesc['layers'], imagedesc['rows'], imagedesc['cols']);
+
+	# the mat (voxel to physical) matrix
 	M = N.eye(4, dtype=N.float64);
-	# for the test data, set the xyz voxel sizes for fMRI volume
-	M[0][0] = 3.75
-	M[1][1] = 3.75
-	M[2][2] = 5.0
+	# for now just the sample size (mm units) in x, y and z
+	M[0][0] = imagedesc['sample_x']
+	M[1][1] = imagedesc['sample_y']
+	M[2][2] = imagedesc['sample_z']
 	# 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
-	max = ImageVolume.max()
-	min = ImageVolume.min()
-	ih  = N.zeros(max-min+1, dtype=N.float64);
-	h   = N.zeros(max-min+1, dtype=N.float64);
-	if threshold <= 0:
-	    threshold = 0.999
-	elif threshold > 1.0:
-	    threshold = 1.0
-	# get the integrated histogram of the volume and get max from 
-	# the threshold crossing in the integrated histogram 
-	index  = R.register_image_threshold(ImageVolume, h, ih, threshold)
-	scale  = 255.0 / (index-min)
-	# generate the scaled 8 bit image
-	images = (scale*(ImageVolume.astype(N.float)-min))
-	images[images>255] = 255 
-	# the data type is now uchar
-	image = {'data' : images.astype(N.uint8), 'mat' : M, 'dim' : D, 'fwhm' : F}
-	if debug == 1:
-    	    return image, h, ih, index
-        else:
+	D[0] = imagedesc['rows']
+	D[1] = imagedesc['cols']
+	D[2] = imagedesc['layers']
+
+	if imagename == None:
+	    # no voxels to scale to 8 bits
+    	    ImageVolume = ImageVolume.astype(N.uint8)
+	    image = {'data' : ImageVolume, 'mat' : M, 'dim' : D, 'fwhm' : F}
     	    return image
 
-
-def load_anatMRI_image(imagename=filename, rows=256, cols=256, layers=90, threshold=0.999, debug=0):
-	# un-scaled images
-    	ImageVolume = N.fromfile(imagename, dtype=N.uint16).reshape(layers, rows, cols);
-	M = N.eye(4, dtype=N.float64);
-	# for the test data, set the xyz voxel sizes for anat-MRI volume
-	M[0][0] = 0.9375
-	M[1][1] = 0.9375
-	M[2][2] = 1.5
-	# 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
+	# 8 bit scale with threshold clip of the volume integrated histogram
 	max = ImageVolume.max()
 	min = ImageVolume.min()
 	ih  = N.zeros(max-min+1, dtype=N.float64);
@@ -481,78 +480,150 @@
 	# generate the scaled 8 bit image
 	images = (scale*(ImageVolume.astype(N.float)-min))
 	images[images>255] = 255 
-	# the data type is now uchar
 	image = {'data' : images.astype(N.uint8), 'mat' : M, 'dim' : D, 'fwhm' : F}
 	if debug == 1:
     	    return image, h, ih, index
         else:
     	    return image
 
+def load_anatMRI_desc():
+	# this is for demo on the test MRI and fMRI volumes
+	rows   = 256
+	cols   = 256
+	layers = 90
+	xsamp  = 0.9375
+	ysamp  = 0.9375
+	zsamp  = 1.5
+	desc = {'rows' : rows, 'cols' : cols, 'layers' : layers, 
+		'sample_x' : xsamp, 'sample_y' : ysamp, 'sample_z' : zsamp}
+	return desc
 
-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);
-	# for the test data, set the xyz voxel sizes for anat-MRI volume
-	M[0][0] = 0.9375
-	M[1][1] = 0.9375
-	M[2][2] = 1.5
-	# 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_fMRI_desc():
+	# this is for demo on the test MRI and fMRI volumes
+	rows   = 64
+	cols   = 64
+	layers = 28
+	xsamp  = 3.75
+	ysamp  = 3.75
+	zsamp  = 5.0
+	desc = {'rows' : rows, 'cols' : cols, 'layers' : layers, 
+		'sample_x' : xsamp, 'sample_y' : ysamp, 'sample_z' : zsamp}
+	return desc
 
 def read_fMRI_directory(path):
 	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
 
-def get_test_rotated_images(alpha=0.0, beta=0.0, gamma=0.0):
-	image1 = load_anatMRI_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
+	ave_fMRI_volume = ave_fMRI_volume / float(count)
 
-def get_test_scaled_images(scale=4.0):
-	# this is for coreg MRI / fMRI test
-	image1 = load_anatMRI_image()
-	image2 = build_scale_image(image1, scale)
-	imdata = build_structs(step=1)
-	# allow the G image to be rotated for testing
-	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
+	return ave_fMRI_volume 
 
+#
+#  ---- testing/debug routines  ---- 
+#
 
 def build_scale_image(image, scale):
 	(layers, rows, cols) = image['data'].shape
-	image2 = image['data'][0:layers:scale, 0:rows:scale, 0:cols:scale]
+	#
+	# rescale the 'mat' (voxel to physical mapping matrix) 
+	#
 	M = image['mat'] * scale
 	# dimensions 
 	D = N.zeros(3, dtype=N.int32);
 	# Gaussian kernel - fill in with build_fwhm() 
 	F = N.zeros(3, dtype=N.float64);
+	Z = N.zeros(3, dtype=N.float64);
 	D[0] = rows/scale
 	D[1] = cols/scale
 	D[2] = layers/scale
+    	image2 = N.zeros(D[2]*D[1]*D[0], dtype=N.uint8).reshape(D[2], D[0], D[1]);
+	mode = 1;
+	R.register_volume_resample(image['data'], image2, Z, scale, mode)
 	scaled_image = {'data' : image2, 'mat' : M, 'dim' : D, 'fwhm' : F}
 	return scaled_image
 
+
+def get_test_MRI_volumes(scale=2, alpha=3.0, beta=4.0, gamma=5.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.  
+	#
+
+	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)
+	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
+	M = build_rotate_matrix(imdata['parms'])
+	R.register_linear_resample(image1['data'], image2['data'], M, imdata['step'])
+	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):
+	#
+	# return rotated fMRIVol
+	#
+
+	desc = load_fMRI_desc()
+	image = load_volume(desc, imagename=None)
+	imdata = build_structs(step=1)
+	image['fwhm'] = build_fwhm(image['mat'], imdata['step'])
+	imdata['parms'][0] = alpha
+	imdata['parms'][1] = beta
+	imdata['parms'][2] = gamma
+	M = build_rotate_matrix(imdata['parms'])
+	R.register_cubic_resample(fMRIVol['data'], image['data'], M, imdata['step'])
+	return image
+
+
+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_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)
+
+	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
+
+



More information about the Scipy-svn mailing list