[Scipy-svn] r3957 - trunk/scipy/ndimage
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Feb 22 20:43:08 EST 2008
Author: tom.waite
Date: 2008-02-22 19:43:06 -0600 (Fri, 22 Feb 2008)
New Revision: 3957
Modified:
trunk/scipy/ndimage/_registration.py
Log:
More fMRI coregistration support added.
Modified: trunk/scipy/ndimage/_registration.py
===================================================================
--- trunk/scipy/ndimage/_registration.py 2008-02-23 01:41:56 UTC (rev 3956)
+++ trunk/scipy/ndimage/_registration.py 2008-02-23 01:43:06 UTC (rev 3957)
@@ -6,6 +6,7 @@
import scipy.ndimage as NDI
import scipy.optimize as OPT
import time
+import glob
# Issue warning regarding heavy development status of this module
import warnings
@@ -72,20 +73,6 @@
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
@@ -209,7 +196,7 @@
Tx=0.0, Ty=0.0, Tz=0.0, stepsize=1):
# takes an image and 3D rotate using trilinear interpolation
- image1 = load_image()
+ image1 = load_anatMRI_image()
image2 = load_blank_image()
imdata = build_structs(step=stepsize)
imdata['parms'][0] = alpha
@@ -227,6 +214,7 @@
image_F_xyz2 = filter_image_3D(image2['data'], image2['fwhm'], ftype)
image2['data'] = image_F_xyz2
+ # this is now a rotated and low pass filtered version of the volume
R.register_linear_resample(image1['data'], image2['data'], M, imdata['step'])
return image2
@@ -246,42 +234,6 @@
# 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]
@@ -461,7 +413,146 @@
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);
+ 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
+ # 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:
+ 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
+ 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:
+ 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);
+ # 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 read_fMRI_directory(path):
+ files_fMRI = glob.glob(path)
+ return files_fMRI
+
+
+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
+
+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
+
+
+def build_scale_image(image, scale):
+ (layers, rows, cols) = image['data'].shape
+ image2 = image['data'][0:layers:scale, 0:rows:scale, 0:cols:scale]
+ 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);
+ D[0] = rows/scale
+ D[1] = cols/scale
+ D[2] = layers/scale
+ scaled_image = {'data' : image2, 'mat' : M, 'dim' : D, 'fwhm' : F}
+ return scaled_image
+
More information about the Scipy-svn
mailing list