[Scipy-svn] r4449 - trunk/scipy/ndimage
scipy-svn@scip...
scipy-svn@scip...
Wed Jun 18 17:36:34 CDT 2008
Author: tom.waite
Date: 2008-06-18 17:36:32 -0500 (Wed, 18 Jun 2008)
New Revision: 4449
Modified:
trunk/scipy/ndimage/_registration.py
Log:
remove demo methods which go to nipy registration
Modified: trunk/scipy/ndimage/_registration.py
===================================================================
--- trunk/scipy/ndimage/_registration.py 2008-06-18 22:15:03 UTC (rev 4448)
+++ trunk/scipy/ndimage/_registration.py 2008-06-18 22:36:32 UTC (rev 4449)
@@ -856,7 +856,7 @@
return rot_matrix
-def build_test_volume(imagedesc, S=[1500.0, 2500.0, 1000.0]):
+def build_gauss_volume(imagedesc, S=[1500.0, 2500.0, 1000.0]):
"""
build a 3D Gaussian volume. user passes in image dims in imagedesc
@@ -900,65 +900,6 @@
return volume3D
-
-def load_volume(imagedesc, imagename=None):
-
- """
-
- returns an nd_array from the filename or blank image. used for testing.
-
- Parameters
- ----------
- imagedesc : {dictionary}
- imagedesc is the descriptor of the image to be read.
-
- imagename : {string} : optional
- name of image file. No name creates a blank image that is used for creating
- a rotated test image or image rescaling.
-
- Returns
- -------
- image : {nd_array}
- the volume data assoicated with the filename or a blank volume of the same
- dimensions as specified in imagedesc.
-
- M : {nd_array}
- the voxel-to-physical affine matrix (mat)
-
- Examples
- --------
-
- >>> import _registration as reg
- >>> anat_desc = reg.load_anatMRI_desc()
- >>> image, M = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img')
-
-
- """
-
- # 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
- image = np.zeros([imagedesc['layers'],imagedesc['rows'],imagedesc['cols']],dtype=np.uint16)
- else:
- image = np.fromfile(imagename,
- dtype=np.uint16).reshape(imagedesc['layers'], imagedesc['rows'], imagedesc['cols']);
-
- # the mat (voxel to physical) matrix
- M = np.eye(4, dtype=np.float64);
- # 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']
-
- if imagename == None:
- # no voxels to scale to 8 bits
- image = image.astype(np.uint8)
-
- return image, M
-
-
-
def scale_image(image, max_amp=255, image_type=np.uint8, threshold=0.999, fetch_ih=0):
"""
@@ -1081,40 +1022,6 @@
-#
-# ---- demo/debug routines ----
-#
-
-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_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_scale_volume(image, mat, scale):
#
# rescale the 'mat' (voxel to physical mapping matrix)
@@ -1133,263 +1040,5 @@
return image2, M
-def demo_build_dual_volumes(scale=2, alpha=3.0, beta=4.0, gamma=5.0, Tx = 0.0, Ty = 0.0, Tz = 0.0):
- """
- demo with (must have file ANAT1_V0001.img)
- builds a volume and a scaled-rotated version for coreg testing
- image1, mat1, image2, mat2 = reg.demo_build_dual_volumes()
- x = reg.register(image1, mat1, image2, mat2, method='ncc', lite=1)
- image2r = reg.remap_image(image2, x, resample='cubic')
- image2rz = reg.resize_image(image2r, mat1)
- """
- #
- # 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.
- #
-
- step = np.array([1, 1, 1], dtype=np.int32)
- anat_desc = load_anatMRI_desc()
- image1, mat1 = load_volume(anat_desc, imagename='ANAT1_V0001.img')
- image2, mat2 = load_volume(anat_desc, imagename=None)
- image1 = scale_image(image1)
- P = np.zeros(6, dtype=np.float64);
- P[0] = alpha
- P[1] = beta
- P[2] = gamma
- P[3] = Tx
- P[4] = Ty
- P[5] = Tz
- M = build_rotate_matrix(P)
- # rotate volume. linear interpolation means the volume is low pass filtered
- reg.register_linear_resample(image1, image2, M, step)
- # subsample volume
- image2, mat2 = build_scale_volume(image2, mat2, scale)
- return image1, mat1, image2, mat2
-
-def tests(image1, mat1, image2, mat2):
-
- # for same image, using the lite method the off-diagonal is zero
- cost, joint = reg.check_alignment(image1, mat1, image2, mat2, ret_histo=1, lite=1)
- my_diag = joint.diagonal()
- Z = np.diag(my_diag)
- W = joint - Z
- W[abs(W) < 1e-10] = 0.0
-
- if W.max() != 0.0 and W.min() != 0.0:
- print 'joint histogram is not diagonal '
- if abs(cost) < 0.99:
- print 'cross correlation is too small'
-
- # for same image, not using the lite method the off-diagonal is non-zero
- cost, joint = reg.check_alignment(image1, mat1, image2, mat2, ret_histo=1, lite=0)
- my_diag = joint.diagonal()
- Z = np.diag(my_diag)
- W = joint - Z
- W[abs(W) < 1e-10] = 0.0
-
- if W.max() == 0.0 and W.min() == 0.0:
- print 'joint histogram is diagonal and needs off-diagonals'
- if abs(cost) < 0.99:
- print 'cross correlation is too small'
-
- # call w/o returning the joint histogram
- cost = reg.check_alignment(image1, mat1, image2, mat2, ret_histo=0, lite=1)
- if abs(cost) < 0.99:
- print 'cross correlation is too small'
-
- cost = reg.check_alignment(image1, mat1, image2, mat2, ret_histo=0, lite=0)
- if abs(cost) < 0.99:
- print 'cross correlation is too small'
-
-
- image1 = np.zeros([64,64,64],np.uint8)
- image2 = np.zeros([64,64,64],np.uint8)
- image3 = np.zeros([64,64],np.uint8)
- mat1 = np.eye(4)
- mat2 = np.eye(3)
- mat3 = np.zeros([4,4])
- # test with wrong dim image, wrong dim mat and mat with zeros on diagonal
-
- # wrong image dim
- assertRaises(ValueError, check_alignment, image1, mat1, image3, mat1)
- # wrong mat dim
- assertRaises(ValueError, check_alignment, image1, mat1, image2, mat2)
- # mat with zeros on diagonal
- assertRaises(ValueError, check_alignment, image1, mat1, image2, mat3)
-
-
-
-
-
-
-
-
-
-
-
-
-def demo_rotate_fMRI_volume(fMRI_volume, desc, x):
- #
- # return rotated fMRIVol.
- #
-
- image = load_volume(desc, imagename=None)
- image = scale_image(image)
- step = np.array([1, 1, 1], dtype=np.int32)
- M = build_rotate_matrix(x)
- # rotate volume. cubic spline interpolation means the volume is NOT low pass filtered
- reg.register_cubic_resample(fMRI_volume, image, M, step)
-
- return image
-
-def demo_MRI_coregistration(anatfile, funclist, optimizer_method='powell',
- histo_method=1, smooth_histo=0, smooth_image=0,
- ftype=1):
- """
- demo with (must have file ANAT1_V0001.img and fMRI directory fMRIData)
-
- measures, imageF_anat, fmri_series = reg.demo_MRI_coregistration()
-
- show results with
-
- In [59]: measures[25]['cost']
- Out[59]: -0.48607185
-
- In [60]: measures[25]['align_cost']
- Out[60]: -0.99514639
-
- In [61]: measures[25]['align_rotate']
- Out[61]:
- array([ 1.94480181, 5.64703989, 5.35002136, -5.00544405, -2.2712214, -1.42249691], dtype=float32)
-
- In [62]: measures[25]['rotate']
- Out[62]:
- array([ 1.36566341, 4.70644331, 4.68198586, -4.32256889, -2.47607017, -2.39173937], dtype=float32)
-
-
- """
-
- # 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.
-
- # read the anatomical MRI volume
- anat_desc = load_anatMRI_desc()
- imageF_anat, anat_mat = load_volume(anat_desc, imagename=anatfile)
- imageF = imageF_anat.copy()
- # the sampling structure
- step = np.array([1, 1, 1], dtype=np.int32)
- # the volume filter
- imageF_anat_fwhm = build_fwhm(anat_mat, step)
-
-
- # allocate the structure for the processed fMRI array
- metric_test = np.dtype([('cost', 'f'),
- ('align_cost', 'f'),
- ('rotate', 'f', 6),
- ('align_rotate', 'f', 6)])
- # allocate the empty dictionary that will contain metrics and aligned volumes
- fmri_series = {}
-
- # read in the file list of the fMRI data
- fMRIdata = read_fMRI_directory('fMRIData\*.img')
- fmri_desc = load_fMRI_desc()
- image_fmri, fmri_mat = load_volume(fmri_desc, fMRIdata[0])
-
- # one time build of the fwhm that is used to build the filter kernels
- anat_fwhm = build_fwhm(anat_mat, step)
- fmri_fwhm = build_fwhm(fmri_mat, step)
-
- # blank volume that will be used for ensemble average for fMRI volumes
- # prior to functional-anatomical coregistration
-<<<<<<< .mine
- ave_fMRI_volume = np.zeros(fmri_desc['layers']*fmri_desc['rows']*fmri_desc['cols'], dtype=np.float64)
-=======
- ave_fMRI_volume = np.zeros([fmri_desc['layers']*fmri_desc['rows']*fmri_desc['cols']],
- dtype=np.float64)
->>>>>>> .r4446
-
- count = 0
- number_volumes = len(fMRIdata)
- measures = np.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:
- fmri_series[count] = image
- count = count + 1
- else:
- x = np.random.random(6) - 0.5
- x = 10.0 * x
- fmri_series[count] = demo_rotate_fMRI_volume(image, fmri_desc, x)
- measures[count]['rotate'][0:6] = x[0:6]
- count = count + 1
-
-
- # load and register the fMRI volumes with volume_0 using normalized cross correlation metric
- imageF = fmri_series[0]
- if smooth_image:
- imageF = filter_image_3D(imageF, fmri_fwhm, ftype)
- for i in range(1, number_volumes):
- imageG = fmri_series[i]
- if smooth_image:
- imageG = filter_image_3D(imageG, fmri_fwhm, ftype)
- # the measure prior to alignment
- measures[i]['cost'] = check_alignment(imageF, fmri_mat, imageG, fmri_mat, method='ncc',
- lite=histo_method, smhist=smooth_histo)
- x = register(imageF, fmri_mat, imageG, fmri_mat, lite=histo_method, method='ncc',
- opt_method=optimizer_method, smhist=smooth_histo)
- measures[i]['align_rotate'][0:6] = x[0:6]
- measures[i]['align_cost'] = check_alignment(imageF, fmri_mat, imageG, fmri_mat,
- 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(np.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(np.float64)
-
- ave_fMRI_volume = (ave_fMRI_volume / float(number_volumes)).astype(np.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:
- imageF_anat = filter_image_3D(imageF_anat, anat_fwhm, ftype)
-
- x = register(imageF_anat, anat_mat, ave_fMRI_volume, fmri_mat, lite=histo_method,
- method='nmi', opt_method=optimizer_method, smhist=smooth_histo)
-
- 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, fmri_series
-
-
-def demo_fMRI_resample(imageF_anat, imageF_anat_mat, 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