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

scipy-svn@scip... scipy-svn@scip...
Mon Jun 9 17:42:05 CDT 2008


Author: tom.waite
Date: 2008-06-09 17:42:02 -0500 (Mon, 09 Jun 2008)
New Revision: 4421

Modified:
   trunk/scipy/ndimage/_registration.py
Log:
bug fixes


Modified: trunk/scipy/ndimage/_registration.py
===================================================================
--- trunk/scipy/ndimage/_registration.py	2008-06-09 06:05:10 UTC (rev 4420)
+++ trunk/scipy/ndimage/_registration.py	2008-06-09 22:42:02 UTC (rev 4421)
@@ -31,69 +31,57 @@
 #  ---- co-registration and IO  ---- 
 #
 
-def resize_image(imageG, imageF_mat):
+def resize_image(imageS, imageS_mat, imageR_mat):
     """
-    zoom_image = resize_image(source_image, reference_image[mat])
+    zoom_image = resize_image(imageS, imageS_mat, imageR_mat)
 
-    Fractional resample source_image to reference_imagesize. The
-    resample is implemented with 3D cubic spline. The reference
-    image [mat] is the 4x4 voxel-to-physical conversion matrix.
+    Fractional resample source_image to reference_image size. The
+    resample is implemented with 3D cubic spline. The source
+    imageS_mat is the 4x4 voxel-to-physical conversion matrix.
     
     Parameters 
     ----------
 
-    imageG : {dictionary} 
-        imageG is the source image to be resized. it is a dictionary with
-        the data as an ndarray in the ['data'] component.
+    imageS: {ndarray} 
+        imageS is the source image to be resized.
 
-    reference_image[mat] : {ndarray}
-        refernce_image is the image whose sampling dimensions the source
-        image is to be remapped to. [mat] refers to the component
-        of the image dictionary, reference_image['mat'] that is the
-        sampling dimensions.
+    imageS_mat : {ndarray} 
+        the 4x4 transform of the source image that maps voxel to physical.
 
+    imageR_mat : {ndarray}
+        the 4x4 transform of the destination image that maps voxel to physical.
+
     Returns 
     -------
-    zoom_image : {dictionary}
+    zoom_image : {ndarray}
 
     Examples
     --------
 
     >>> import _registration as reg
-    >>> measures, imageF_anat, fmri_series = reg.demo_MRI_coregistration()
+    >>> measures, image_anat, image_anat_mat, image_fmri_mat, fmri_series = reg.demo_MRI_coregistration()
 
-    >>> resampled_fmri = reg.resize_image(fmri_series[10], imageF_anat['mat'])
+    >>> resampled_fmri = reg.resize_image(fmri_series[10], image_fmri_mat, image_anat_mat)
 
-    image 10 in the fmri_series is resampled to imageF_anat coordinates
+    image 10 in the fmri_series is resampled from image_fmri_mat to image_anat coordinates
 
     """
 
-    Z = np.zeros(3, dtype=np.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]
+    Z = imageS.diagonal() / imageR.diagonal()
 
-    # new volume dimensions (rounded)
-    D = np.zeros(3, dtype=np.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)
+    # new volume dimensions (rounded). D, imageS and Z are 3D and this is a vector element product
+    D = (imageS.shape * Z + 0.5).astype(np.int16)
 
-    M = np.eye(4, dtype=np.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 = np.zeros(D[2]*D[1]*D[0], dtype=np.uint8).reshape(D[2], D[0], D[1])
+    # for the test data, set the xyz voxel sizes for fMRI volume. M is a 4x4 matrix.
+    M = np.diag(imageS.diagonal() / Z)    
+    image = np.empty((D[2],D[1],D[0]),np.uint8)
+    
     mode  = 2
     scale = 0
-    reg.register_volume_resample(imageG['data'], image, Z, scale, mode)
-    F = np.zeros(3, dtype=np.float64);
-    zoom_image = {'data' : image, 'mat' : M, 'dim' : D, 'fwhm' : F}
+    reg.register_volume_resample(imageS, image, Z, scale, mode)
 
-    return zoom_image
+    return image, M
 
 def remap_image(image, parm_vector, resample='linear'):
     """
@@ -105,20 +93,19 @@
 
     Parameters 
     ----------
-    image : {dictionary} 
-        image is the source image to be remapped. it is a dictionary with
-        the data as an ndarray in the ['data'] component.
+    image : {ndarray} 
+        image is the source image to be remapped. 
 
     parm_vector : {ndarray}
         parm_vector is the 6-dimensional vector (3 angles, 3 translations)
-        generated from the registration.
+        generated from the rigid body registration.
 
     resample : {'linear', 'cubic'}, optional
 
 
     Returns 
     -------
-    remaped_image : {dictionary}
+    remaped_image : {ndarray}
 
     Examples
     --------
@@ -134,19 +121,20 @@
     # 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
+    (layers, rows, cols) = image.shape
+
     # allocate the zero image
-    remaped_image = np.zeros(layers*rows*cols, dtype=np.uint8).reshape(layers, rows, cols)
-    remaped_image = {'data' : remaped_image, 'mat' : image['mat'], 
-                     'dim' : image['dim'], 'fwhm' : image['fwhm']}
-    imdata = build_structs()
+    #remaped_image = np.zeros(image.size, dtype=np.uint8).reshape(layers, rows, cols)
+    remaped_image = np.empty(image.shape, dtype=np.uint8)
 
+    step = np.array([1, 1, 1], dtype=np.int32)
+
     if resample == 'linear':
         # trilinear interpolation mapping.
-        reg.register_linear_resample(image['data'], remaped_image['data'], M_inverse, imdata['step'])
+        reg.register_linear_resample(image, remaped_image, M_inverse, step)
     elif resample == 'cubic':
         # tricubic convolve interpolation mapping. 
-        reg.register_cubic_resample(image['data'], remaped_image['data'], M_inverse, imdata['step'])
+        reg.register_cubic_resample(image, remaped_image, M_inverse, step)
 
     return remaped_image
 
@@ -182,24 +170,18 @@
 
     """
     # get the inverse mapping to rotate the G matrix to F space following registration
-    imdata = build_structs()
-    # inverse angles and translations
-    imdata['parms'][0] = -parm_vector[0]
-    imdata['parms'][1] = -parm_vector[1]
-    imdata['parms'][2] = -parm_vector[2]
-    imdata['parms'][3] = -parm_vector[3]
-    imdata['parms'][4] = -parm_vector[4]
-    imdata['parms'][5] = -parm_vector[5]
-    M_inverse = build_rotate_matrix(imdata['parms'])
+    # -parm_vector is the inverse angles and translations
+    M_inverse = build_rotate_matrix(-parm_vector)
     return M_inverse
 
-def python_coreg(image1, image2, imdata, ftype=1, smimage=0, lite=0, smhist=0,
-                 method='nmi', opt_method='powell'):
+def coregister(image1, image1_mat, image2, image2_mat, multires=[4, 2], histo_fwhm=3, 
+               ftype=1, lite=0, smhist=0, method='nmi', opt_method='powell'):
+
     """
-    parm_vector = python_coreg(image1, image2, imdata, ftype=1, smimage=0, lite=0,
-                               smhist=0, method='nmi', opt_method='powell'):
+    parm_vector = coregister(image1, image1_mat, image2, image2_mat, multires=[4, 2], histo_fwhm=3,
+                             ftype=1, lite=0, smhist=0, method='nmi', opt_method='powell'):
 
-    takes two images and the image data descriptor (imdata) and determines the optimal 
+    takes two images and the image process descriptor (improc) and determines the optimal 
     alignment of the two images (measured by mutual information or cross correlation) 
     using optimization search of 3 angle and 3 translation parameters. The optimization 
     uses either the Powell or Conjugate Gradient methods in the scipy optimization 
@@ -207,29 +189,33 @@
 
     Parameters 
     ----------
-    image1 : {dictionary} 
+    image1 : {nd_array} 
         image1 is the source image to be remapped during the registration. 
-        it is a dictionary with the data as an ndarray in the ['data'] component.
-    image2 : {dictionary} 
+    image1_mat : {nd_array} 
+        image1_mat is the source image MAT 
+    image2 : {nd_array} 
         image2 is the reference image that image1 gets mapped to. 
-    imdata : {dictionary} 
-        image sampling and optimization information.
+    image2_mat : {nd_array} 
+        image2_mat is the source image MAT 
+    multires: {list}, optional
+        the volume subsample values for each pass of the registration.
+	the default is 2 passes with subsample 4 in pass 1 and subsample 2 in pass 2
+    histo_fwhm : {int}, optional
+        used for the filter kernel in the low pass filter of the joint histogram 
     ftype : {0, 1}, optional
         flag for type of low pass filter. 0 is Gauss-Spline
         1 is pure Gauss. Sigma determined from volume sampling info.
-    smimage : {0, 1}, optional
-        flag for volume 3D low pass filtering of image 2.
-        0 for no filter, 1 for do filter.
     lite : {0, 1}, optional
         lite of 1 is to jitter both images during resampling. 0
         is to not jitter. jittering is for non-aliased volumes.
     smhist: {0, 1}, optional
         flag for joint histogram low pass filtering. 0 for no filter,
         1 for do filter.
-    method: {'nmi', 'mi', 'ncc', 'ecc'}, optional
+    method: {'nmi', 'mi', 'ncc', 'ecc', 'mse'}, optional
         flag for type of registration metric. nmi is normalized mutual
         information; mi is mutual information; ecc is entropy cross
-        correlation; ncc is normalized cross correlation.
+        correlation; ncc is normalized cross correlation. mse is mean
+	squared error.
     opt_method: {'powell', 'hybrid'}, optional
         registration is two pass. Pass 1 is low res to get close to alignment
         and pass 2 starts at the pass 1 optimal alignment. In powell pass 1 and
@@ -248,34 +234,33 @@
     >>> import numpy as NP
     >>> import _registration as reg
 
-    >>> image1, image2, imdata = reg.demo_MRI_volume_align()
-    >>> parm_vector = python_coreg(image1, image2, imdata)
+    >>> image1, image2, fwhm, improc = reg.demo_build_dual_volumes()
+    >>> parm_vector = coregister(image1, image2, fwhm, improc)
 
     """
+
     start = time.time()
-    # smooth of the images
-    if smimage: 
-        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)
+    parm_vector = multires_registration(image1, image1_mat, image2, image2_mat, multires, 
+		                        histo_fwhm, lite, smhist, method, opt_method)
     stop = time.time()
     print 'Total Optimizer Time is ', (stop-start)
     return parm_vector
 
-def multires_registration(image1, image2, imdata, lite, smhist, method, opt_method):
+def multires_registration(image1, image1_mat, image2, image2_mat, multires, histo_fwhm, 
+		          lite, smhist, method, opt_method):
+
     """
     x = multires_registration(image1, image2, imdata, lite, smhist, method, opt_method)
 
-    to be called by python_coreg() which optionally does 3D image filtering and 
+    to be called by coregister() which optionally does 3D image filtering and 
     provies timing for registration.
 
     Parameters 
     ----------
 
-    image1 : {dictionary} 
+    image1 : {nd_array} 
         image1 is the source image to be remapped during the registration. 
-        it is a dictionary with the data as an ndarray in the ['data'] component.
-    image2 : {dictionary} 
+    image2 : {nd_array} 
         image2 is the reference image that image1 gets mapped to. 
     imdata : {dictionary} 
         image sampling and optimization information.
@@ -285,10 +270,11 @@
     smhist: {integer}
         flag for joint histogram low pass filtering. 0 for no filter,
         1 for do filter.
-    method: {'nmi', 'mi', 'ncc', 'ecc'}
+    method: {'nmi', 'mi', 'ncc', 'ecc', 'mse'}
         flag for type of registration metric. nmi is normalized mutual
         information; mi is mutual information; ecc is entropy cross
-        correlation; ncc is normalized cross correlation.
+        correlation; ncc is normalized cross correlation. mse is mean
+	square error.
     opt_method: {'powell', 'hybrid'}
         registration is two pass. Pass 1 is low res to get close to alignment
         and pass 2 starts at the pass 1 optimal alignment. In powell pass 1 and
@@ -303,31 +289,32 @@
     Examples
     --------
 
-    (calling this from python_coreg which optionally filters image2)
+    (calling this from coregister which optionally filters image2)
     >>> import numpy as NP
     >>> import _registration as reg
-    >>> image1, image2, imdata = reg.demo_MRI_volume_align()
-    >>> parm_vector = python_coreg(image1, image2, imdata)
+    >>> image1, mat1, image2, mat2 = reg.demo_build_dual_volumes()
+    >>> parm_vector = coregister(image1, image2, imdata)
 
     """
     ret_histo=0
-    # zero out the start parameter; but this may be set to large values 
-    # if the head is out of range and well off the optimal alignment skirt
-    imdata['parms'][0:5] = 0.0
+    step = np.array([1, 1, 1], dtype=np.int32)
+    fwhm = np.zeros(2, dtype=np.int32)
     # make the step a scalar to can put in a multi-res loop
-    loop = range(imdata['sample'].size)
-    x = imdata['parms']
+    loop = range(size(multires))
+    # 6-D zero vector
+    x = np.zeros(6, dtype=np.float64);
+    # the kernel fwhm value for the x and y joint histogram filter
+    fwhm[:] = histo_fwhm
     for i in loop:
-        step = imdata['sample'][i]
-        imdata['step'][:] = step
-        optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist,
-                        method, ret_histo)
+	# this is the volume subsample
+	step[:] = multires[i]
+        optfunc_args = (image1, image1_mat, image2, image2_mat, step, fwhm, lite,
+			smhist, method, ret_histo)
         p_args = (optfunc_args,)
         if opt_method=='powell':
             print 'POWELL multi-res registration step size ', step
             print 'vector ', x
-            x = fmin_powell(optimize_function, x, args=p_args,
-                                callback=callback_powell) 
+            x = fmin_powell(optimize_function, x, args=p_args, callback=callback_powell)
         elif opt_method=='cg':
             print 'CG multi-res registration step size ', step
             print 'vector ', x
@@ -337,16 +324,16 @@
                 print 'Hybrid POWELL multi-res registration step size ', step
                 print 'vector ', x
                 lite = 0
-                optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist,
-                                method, ret_histo)
+                optfunc_args = (image1, image1_mat, image2, image2_mat, step, fwhm, lite,
+                                smhist, method, ret_histo)
                 p_args = (optfunc_args,)
                 x = fmin_powell(optimize_function, x, args=p_args, callback=callback_powell) 
             elif i==1:
                 print 'Hybrid CG multi-res registration step size ', step
                 print 'vector ', x
                 lite = 1
-                optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, 
-                                smhist, method, ret_histo)
+                optfunc_args = (image1, image1_mat, image2, image2_mat, step, histo_fwhm,
+                                lite, smhist, method, ret_histo)
                 p_args = (optfunc_args,)
                 x = fmin_cg(optimize_function, x, args=p_args, callback=callback_cg) 
 
@@ -369,7 +356,7 @@
     print x
     return
 
-def smooth_kernel(fwhm, x, ktype=1):
+def smooth_kernel(fwhm, x, pixel_scale=8.0, ktype=1):
     """
     kernel = smooth_kernel(fwhm, x, ktype=1)
 
@@ -411,7 +398,7 @@
 
     """
     eps = 0.00001
-    s   = np.square((fwhm/math.sqrt(8.0*math.log(2.0)))) + eps
+    s   = np.square((fwhm/math.sqrt(pixel_scale*math.log(2.0)))) + eps
     if ktype==1:
         # from SPM: Gauss kernel convolved with 1st degree B spline
         w1 = 0.5 * math.sqrt(2.0/s)
@@ -453,7 +440,7 @@
     --------
 
     >>> import _registration as reg
-    >>> image1, image2, imdata = reg.demo_MRI_volume_align()
+    >>> image1, image2, imdata = reg.demo_build_dual_volumes()
     >>> ftype = 1
     >>> image_Filter_xyz = filter_image_3D(image1['data'], image1['fwhm'], ftype)
     >>> image1['data'] = image_Filter_xyz
@@ -462,12 +449,15 @@
     p = np.ceil(2*fwhm[0]).astype(int)
     x = np.array(range(-p, p+1))
     kernel_x = smooth_kernel(fwhm[0], x, ktype=ftype)
+
     p = np.ceil(2*fwhm[1]).astype(int)
     x = np.array(range(-p, p+1))
     kernel_y = smooth_kernel(fwhm[1], x, ktype=ftype)
+
     p = np.ceil(2*fwhm[2]).astype(int)
     x = np.array(range(-p, p+1))
     kernel_z = smooth_kernel(fwhm[2], x, ktype=ftype)
+
     output=None
     # 3D filter in 3 1D separable stages
     axis = 0
@@ -504,14 +494,14 @@
     >>> import _registration as reg
     >>> anat_desc = reg.load_anatMRI_desc()
     >>> image1 = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img')
-    >>> imdata = reg.build_structs()
     >>> image1['fwhm'] = reg.build_fwhm(image1['mat'], imdata['step'])
 
     """
+    # M contains the voxel to physical mapping
     view_3x3 = np.square(M[0:3, 0:3])
-    # sum the elements inn the first row
+    # sum the elements in the first row
     vxg = np.sqrt(view_3x3.sum(axis=0))
-    # assumes that sampling is the same for xyz
+    # assumes that voxel sampling is the same for xyz as S is the step
     size = np.array([1,1,1])*S[0]
     x = np.square(size) - np.square(vxg)
     # clip
@@ -592,7 +582,6 @@
     >>> anat_desc = reg.load_anatMRI_desc()
     >>> image1 = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img')
     >>> image2 = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img')
-    >>> imdata = reg.build_structs()
     >>> image1['fwhm'] = reg.build_fwhm(image1['mat'], imdata['step'])
     >>> image2['fwhm'] = reg.build_fwhm(image2['mat'], imdata['step'])
     >>> method = 'ncc'
@@ -607,24 +596,26 @@
     """
 
     image_F       = optfunc_args[0]
-    image_G       = optfunc_args[1]
-    sample_vector = optfunc_args[2]
-    fwhm          = optfunc_args[3]
-    do_lite       = optfunc_args[4]
-    smooth        = optfunc_args[5]
-    method        = optfunc_args[6]
-    ret_histo     = optfunc_args[7]
+    image_F_mat   = optfunc_args[1]
+    image_G       = optfunc_args[2]
+    image_G_mat   = optfunc_args[3]
+    sample_vector = optfunc_args[4]
+    fwhm          = optfunc_args[5]
+    do_lite       = optfunc_args[6]
+    smooth        = optfunc_args[7]
+    method        = optfunc_args[8]
+    ret_histo     = optfunc_args[9]
 
     rot_matrix = build_rotate_matrix(x)
     cost = 0.0
     epsilon = 2.2e-16 
     # image_G is base image
     # image_F is the to-be-rotated image
-    # rot_matrix is the 4x4 constructed (current angles and translates) transform matrix
+    # rot_matrix is the 4x4 constructed (rigid body) transform matrix
     # sample_vector is the subsample vector for x-y-z
 
-    F_inv = np.linalg.inv(image_F['mat'])
-    composite = np.dot(F_inv, image_G['mat'])
+    F_inv = np.linalg.inv(image_F_mat)
+    composite = np.dot(F_inv, image_G_mat)
     composite = np.dot(composite, rot_matrix)
 
     if method == 'mse':
@@ -632,14 +623,14 @@
         # mean squard error method
         #
 
-        (layers, rows, cols) = image_F['data'].shape
         # allocate the zero image
-        remap_image_F = np.zeros(layers*rows*cols, dtype=np.uint8).reshape(layers, rows, cols)
-        imdata = build_structs()
+        #(layers, rows, cols) = image_F.shape
+        #remap_image_F = np.zeros(image_F.size, dtype=np.uint8).reshape(layers, rows, cols)
+        remap_image_F = np.empty(image_F.shape, dtype=np.uint8)
         # trilinear interpolation mapping.
-        reg.register_linear_resample(image_F['data'], remap_image_F, composite,
-                                   imdata['step'])
-        cost = (np.square(image_G['data']-remap_image_F)).mean()
+        reg.register_linear_resample(image_F, remap_image_F, composite, sample_vector)
+        cost = (np.square(image_G-remap_image_F)).mean()
+	# cost is min when G and F are aligned so keep cost positive
 
         return cost
 
@@ -649,29 +640,25 @@
         #
 
         # allocate memory for 2D histogram
-        joint_histogram = np.zeros([256, 256], dtype=np.float64);
+        joint_histogram = np.zeros([256, 256], dtype=np.float64)
 
         if do_lite: 
-            reg.register_histogram_lite(image_F['data'], image_G['data'], composite,
-                                      sample_vector, joint_histogram)
+            reg.register_histogram_lite(image_F, image_G, composite, sample_vector, joint_histogram)
         else:
-            reg.register_histogram(image_F['data'], image_G['data'], composite,
-                                 sample_vector, joint_histogram)
+            reg.register_histogram(image_F, image_G, composite, sample_vector, joint_histogram)
 
         # smooth the histogram
         if smooth: 
-            p = np.ceil(2*fwhm[0]).astype(int)
+            p = np.ceil(2*fwhm).astype(int)
             x = np.array(range(-p, p+1))
-            kernel1 = smooth_kernel(fwhm[0], x)
-            p = np.ceil(2*fwhm[1]).astype(int)
-            x = np.array(range(-p, p+1))
-            kernel2 = smooth_kernel(fwhm[1], x)
+            hkernel = smooth_kernel(fwhm, x)
             output=None
-            # 2D filter in 1D separable stages
+            # 2D filter in 1D separable stages using the same kernel. SPM
+	    # has options for a 2D fwhm kernel yet only uses 1 element
             axis = 0
-            result = correlate1d(joint_histogram, kernel1, axis, output)
+            joint_histogram = correlate1d(joint_histogram, hkernel, axis, output)
             axis = 1
-            joint_histogram = correlate1d(result, kernel1, axis, output)
+            joint_histogram = correlate1d(joint_histogram, hkernel, axis, output)
 
         joint_histogram += epsilon # prevent log(0) 
         # normalize the joint histogram
@@ -702,7 +689,7 @@
             row_entropy = marginal_row * np.log(marginal_row)
             col_entropy = marginal_col * np.log(marginal_col)
             H = joint_histogram * np.log(joint_histogram)  
-            nmi = (row_entropy.sum() + col_entropy.sum()) / (H.sum())
+            nmi  = (row_entropy.sum() + col_entropy.sum()) / (H.sum())
             cost = -nmi
 
         elif method == 'ncc':
@@ -719,7 +706,7 @@
             b = b - m2
             # element multiplies in the joint histogram and grids
             H = ((joint_histogram * a) * b).sum()
-            ncc = H / (np.dot(sig1, sig2)) 
+            ncc  = H / (np.dot(sig1, sig2)) 
             cost = -ncc
 
         if ret_histo:
@@ -728,64 +715,6 @@
             return cost
 
 
-def build_structs(step=1):
-    """
-    img_data = build_structs(step=1)
-
-    builds the image data (imdata) dictionary for later use as parameter
-    storage in the co-registration.
-
-    Parameters 
-    ----------
-    step : {int} : optional
-    default is 1 and is the sample increment in voxels. This sets the sample
-    for x,y,z and is the same value in all 3 axes. only change the default for debug.
-
-    Returns 
-    -------
-    img_data : {dictionary}
-
-    Examples
-    --------
-
-    >>> import numpy as NP
-    >>> import _registration as reg
-    >>> imdata = reg.build_structs()
-
-    """
-
-    # build image data structures here
-    P = np.zeros(6, dtype=np.float64);
-    T = np.zeros(6, dtype=np.float64);
-    F = np.zeros(2, dtype=np.int32);
-    S = np.ones(3,  dtype=np.int32);
-    sample = np.zeros(2, dtype=np.int32);
-    S[0] = step
-    S[1] = step
-    S[2] = step
-    # image/histogram smoothing
-    F[0] = 3
-    F[1] = 3
-    # subsample for multiresolution registration
-    sample[0] = 4
-    sample[1] = 2
-    # tolerances for angle (0-2) and translation (3-5)
-    T[0] = 0.02 
-    T[1] = 0.02 
-    T[2] = 0.02 
-    T[3] = 0.001 
-    T[4] = 0.001 
-    T[5] = 0.001 
-    # P[0] = alpha <=> pitch. + alpha is moving back in the sagittal plane
-    # P[1] = beta  <=> roll.  + beta  is moving right in the coronal plane
-    # P[2] = gamma <=> yaw.   + gamma is right turn in the transverse plane
-    # P[3] = Tx
-    # P[4] = Ty
-    # P[5] = Tz
-    img_data = {'parms' : P, 'step' : S, 'fwhm' : F, 'tol' : T, 'sample' : sample}
-    return img_data
-
-
 def build_rotate_matrix(img_data_parms):
     """
     rot_matrix = reg.build_rotate_matrix(img_data_parms)
@@ -807,7 +736,6 @@
 
     >>> import numpy as NP
     >>> import _registration as reg
-    >>> imdata = reg.build_structs()
     >>> x = np.zeros(6, dtype=np.float64)
     >>> M = reg.build_rotate_matrix(x)
     >>> M 
@@ -864,6 +792,51 @@
     return rot_matrix
 
 
+def build_test_volume(imagedesc, S=[15.0, 25.0, 10.0]):
+
+    """
+    build a 3D Gaussian volume. user passes in image dims in imagedesc
+    the sigma for each axis is S[3] where 0=z, 1=y, 2=x
+
+    volume3D = build_test_volume(imagedesc, S)
+
+    Parameters 
+    ----------
+    imagedesc : {dictionary}
+        volume dimensions and sampling
+
+    S : {tuple}
+        the Gaussian sigma for Z, Y and X
+
+    Returns 
+    -------
+
+    volume3D : {nd_array}
+        the 3D volume for testing
+
+    """
+    layers = imagedesc['layers']
+    rows   = imagedesc['rows']
+    cols   = imagedesc['cols']
+
+    L = layers/2
+    R = rows/2
+    C = cols/2
+
+    # build coordinates for 3D Gaussian volume
+    # coordinates are centered at (0, 0, 0)
+    [a, b, c] = np.mgrid[-L:L, -R:R, -C:C]
+
+    sigma    = np.array([S[0], S[1], S[2]])
+    aa       = (np.square(a))/sigma[0]
+    bb       = (np.square(b))/sigma[1]
+    cc       = (np.square(c))/sigma[2]
+    volume3D = np.exp(-(aa + bb + cc))
+
+    return volume3D
+
+
+
 def load_volume(imagedesc, imagename=None, threshold=0.999, debug=0):
 
     """
@@ -895,13 +868,13 @@
 
     Returns 
     -------
-    image : {dictionary}
+    image : {nd_array}
         the volume data assoicated with the filename or a blank volume of the same
         dimensions as specified in imagedesc.
 
     --- OR --- (if debug = 1)
 
-    image : {dictionary}
+    image : {nd_array}
         the volume data assoicated with the filename or a blank volume of the same
         dimensions as specified in imagedesc.
 
@@ -932,11 +905,12 @@
     # 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 = np.zeros(imagedesc['layers']*imagedesc['rows']*imagedesc['cols'],
-                        dtype=np.uint16).reshape(imagedesc['layers'], imagedesc['rows'], imagedesc['cols'])
+        image = np.zeros([imagedesc['layers'],imagedesc['rows'],imagedesc['cols']],dtype=np.uint16)
+        #image = np.zeros(imagedesc['layers']*imagedesc['rows']*imagedesc['cols'],
+        #           dtype=np.uint16).reshape(imagedesc['layers'], imagedesc['rows'], imagedesc['cols'])
     else:
-        ImageVolume = np.fromfile(imagename,
-                        dtype=np.uint16).reshape(imagedesc['layers'], imagedesc['rows'], imagedesc['cols']);
+        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);
@@ -944,23 +918,15 @@
     M[0][0] = imagedesc['sample_x']
     M[1][1] = imagedesc['sample_y']
     M[2][2] = imagedesc['sample_z']
-    # dimensions 
-    D = np.zeros(3, dtype=np.int32);
-    # Gaussian kernel - fill in with build_fwhm() 
-    F = np.zeros(3, dtype=np.float64);
-    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(np.uint8)
-        image = {'data' : ImageVolume, 'mat' : M, 'dim' : D, 'fwhm' : F}
-        return image
+        image = image.astype(np.uint8)
+        return image, M
 
     # 8 bit scale with threshold clip of the volume integrated histogram
-    max = ImageVolume.max()
-    min = ImageVolume.min()
+    max = image.max()
+    min = image.min()
     ih  = np.zeros(max-min+1, dtype=np.float64);
     h   = np.zeros(max-min+1, dtype=np.float64);
     if threshold <= 0:
@@ -969,16 +935,16 @@
         threshold = 1.0
     # get the integrated histogram of the volume and get max from 
     # the threshold crossing in the integrated histogram 
-    index  = reg.register_image_threshold(ImageVolume, h, ih, threshold)
+    index  = reg.register_image_threshold(image, h, ih, threshold)
     scale  = 255.0 / (index-min)
     # generate the scaled 8 bit image
-    images = (scale*(ImageVolume.astype(np.float)-min))
-    images[images>255] = 255 
-    image = {'data' : images.astype(np.uint8), 'mat' : M, 'dim' : D, 'fwhm' : F}
+    image = (scale*(image.astype(np.float)-min))
+    image[image>255] = 255 
+    image = image.astype(np.uint8)
     if debug == 1:
-        return image, h, ih, index
+        return image, M, h, ih, index
     else:
-        return image
+        return image, M
 
 
 
@@ -1015,120 +981,100 @@
     return files_fMRI
 
 
-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):
+def check_alignment(image1, image1_mat, image2, image2_mat, histo_fwhm=3, 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
+    # to test the cost function and (optional) view the joint histogram
+    # default of use of ncc for testing the cross-correlation as a metric
+    # of alignment
     #
-    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)
 
+    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
+
+    step = np.array([1, 1, 1], dtype=np.int32)
+    optfunc_args = (image1, image1_mat, image2, image2_mat, step, histo_fwhm, lite,
+		    smhist, method, ret_histo)
+			
+    #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)
+        cost, joint_histogram = optimize_function(P, optfunc_args)
         return cost, joint_histogram 
     else:
-        cost = optimize_function(imdata['parms'], optfunc_args)
+        cost = optimize_function(P, optfunc_args)
         return cost
 
-def build_scale_image(image, scale):
+def build_scale_volume(image, mat, scale):
     #
     # rescale the 'mat' (voxel to physical mapping matrix) 
     #
-    (layers, rows, cols) = image['data'].shape
-    M = image['mat'] * scale
+    (layers, rows, cols) = image.shape
+    M = mat * scale
     # dimensions 
     D = np.zeros(3, dtype=np.int32);
-    # Gaussian kernel - fill in with build_fwhm() 
-    F = np.zeros(3, dtype=np.float64);
     Z = np.zeros(3, dtype=np.float64);
     D[0] = rows/scale
     D[1] = cols/scale
     D[2] = layers/scale
-    image2 = np.zeros(D[2]*D[1]*D[0], dtype=np.uint8).reshape(D[2], D[0], D[1]);
+    #image2 = np.zeros(D.prod(), dtype=np.uint8).reshape(D[2], D[0], D[1]);
+    image2 = np.empty([D[2], D[0], D[1]], dtype=np.uint8)
     mode = 1;
-    reg.register_volume_resample(image['data'], image2, Z, scale, mode)
-    scaled_image = {'data' : image2, 'mat' : M, 'dim' : D, 'fwhm' : F}
-    return scaled_image
+    reg.register_volume_resample(image, image2, Z, scale, mode)
+    return image2, M
 
 
-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):
+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, image2, imdata = reg.demo_MRI_volume_align()
-    x = reg.python_coreg(image1, image2, imdata, method='ncc', lite=1) 
+    image1, mat1, image2, mat2 = reg.demo_build_dual_volumes()
+    x = reg.coregister(image1, mat1, image2, mat2, method='ncc', lite=1) 
     image2r = reg.remap_image(image2, x, resample='cubic')
-    image2rz = reg.resize_image(image2r, image1['mat'])
+    image2rz = reg.resize_image(image2r, mat1)
 
-
-    slice1 = image1['data'][45, :, :]
-    slice2 = image2['data'][45/2, :, :]
-    slice2r = image2r['data'][45/2, :, :]
-    slice2rz = image2rz['data'][45, :, :]
-
-    pylab.figure(1)
-    pylab.bone()
-    pylab.imshow(slice1)
-    pylab.imshow(slice1)
-    pylab.figure(2)
-    pylab.imshow(slice2)
-    pylab.figure(3)
-    pylab.imshow(slice2r)
-    pylab.figure(4)
-    pylab.imshow(slice2rz)
-    pylab.show()
-
     """
     #
     # 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 = load_volume(anat_desc, imagename='ANAT1_V0001.img')
-    image2 = load_volume(anat_desc, imagename=None)
-    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'])
+    image1, mat1 = load_volume(anat_desc, imagename='ANAT1_V0001.img')
+    image2, mat2 = load_volume(anat_desc, imagename=None)
+    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['data'], image2['data'], M, imdata['step'])
+    reg.register_linear_resample(image1, image2, M, step)
     # subsample volume
-    image3 = build_scale_image(image2, scale)
-    return image1, image3, imdata
+    image2, mat2 = build_scale_volume(image2, mat2, scale)
+    return image1, mat1, image2, mat2
 
-def demo_rotate_fMRI_volume(fMRIVol, x): 
+def demo_rotate_fMRI_volume(fMRI_volume, desc, x): 
     #
-    # return rotated fMRIVol. the fMRIVol is already loaded, and gets rotated
+    # return rotated fMRIVol.
     #
 
-    desc = load_fMRI_desc()
     image = load_volume(desc, imagename=None)
-    imdata = build_structs()
-    image['fwhm'] = build_fwhm(image['mat'], imdata['step'])
-    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'])
+    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(fMRIVol['data'], image['data'], M, imdata['step'])
+    reg.register_cubic_resample(fMRI_volume, image, M, step)
+
     return image
 
 def demo_MRI_coregistration(anatfile, funclist, optimizer_method='powell', 
@@ -1169,25 +1115,35 @@
 
     # read the anatomical MRI volume
     anat_desc = load_anatMRI_desc()
-    imageF_anat = load_volume(anat_desc, imagename=anatfile)
+    imageF_anat, anat_mat = load_volume(anat_desc, imagename=anatfile)
+    imageF = imageF_anat.copy()
     # the sampling structure
-    imdata = build_structs()
+    step = np.array([1, 1, 1], dtype=np.int32)
     # the volume filter
-    imageF_anat['fwhm'] = build_fwhm(imageF_anat['mat'], imdata['step'])
+    imageF_anat_fwhm = build_fwhm(mat_anat, step)
 
-    # read in the file list of the fMRI data
+
+    # allocate the structure for the processed fMRI array
     metric_test = np.dtype([('cost', 'f'),
-                           ('align_cost', 'f'),
-                           ('rotate', 'f', 6),
-                           ('align_rotate', 'f', 6)])
+                            ('align_cost', 'f'),
+                            ('rotate', 'f', 6),
+                            ('align_rotate', 'f', 6)])
+    # allocate the empty dictionary that will contain metrics and aligned volumes 
+    fmri_series = {}
 
-    #fMRIdata = read_fMRI_directory('fMRIData\*.img')
-    #fMRIdata = read_fMRI_directory(funcdir + '/*.img')
-    fMRIdata = funclist
+    # read in the file list of the fMRI data
+    fMRIdata = read_fMRI_directory('fMRIData\*.img')
     fmri_desc = load_fMRI_desc()
-    fmri_series = {}
-    ave_fMRI_volume = np.zeros(fmri_desc['layers']*fmri_desc['rows']*fmri_desc['cols'],
-                      dtype=np.float64).reshape(fmri_desc['layers'], fmri_desc['rows'], fmri_desc['cols'])
+    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
+    ave_fMRI_volume = np.zeros([fmri_desc['layers']*fmri_desc['rows']*fmri_desc['cols']].dtype=np.float64)
+
     count = 0
     number_volumes = len(fMRIdata)
     measures = np.zeros(number_volumes, dtype=metric_test)
@@ -1196,13 +1152,12 @@
         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 = np.random.random(6) - 0.5
             x = 10.0 * x
-            fmri_series[count] = demo_rotate_fMRI_volume(image, x)
+            fmri_series[count] = demo_rotate_fMRI_volume(image, fmri_desc, x)
             measures[count]['rotate'][0:6] = x[0:6]
             count = count + 1
 
@@ -1210,19 +1165,22 @@
     # 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
+        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, imageG, imdata, method='ncc',
+        measures[i]['cost'] = check_alignment(imageF, fmri_mat, imageG, fmri_mat, 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)
+        x = coregister(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, 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])
+        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 
@@ -1239,10 +1197,11 @@
                        '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)
+        imageF_anat = filter_image_3D(imageF_anat, anat_fwhm, ftype)
+
+    x = coregister(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):
@@ -1250,14 +1209,14 @@
         # overwrite the fMRI volume with the anatomical-aligned volume
         fmri_series[i] = remap_image(image, x, resample='cubic')
 
-    return measures, imageF_anat, fmri_series
+    return measures, imageF, fmri_series
 
 
-def demo_fMRI_resample(imageF_anat, 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'])
+        resampled_fmri_series[i] = resize_image(fmri_series[i], imageF_anat_mat)
 
     return resampled_fmri_series
 



More information about the Scipy-svn mailing list