[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