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

scipy-svn@scip... scipy-svn@scip...
Wed Jun 18 14:04:18 CDT 2008


Author: tom.waite
Date: 2008-06-18 14:04:15 -0500 (Wed, 18 Jun 2008)
New Revision: 4447

Modified:
   trunk/scipy/ndimage/_registration.py
Log:
Parameter checking. Replace c-extension integrated histogram thresholding with pure numpy version.

Modified: trunk/scipy/ndimage/_registration.py
===================================================================
--- trunk/scipy/ndimage/_registration.py	2008-06-17 01:19:04 UTC (rev 4446)
+++ trunk/scipy/ndimage/_registration.py	2008-06-18 19:04:15 UTC (rev 4447)
@@ -74,8 +74,9 @@
     D = (imageS.shape * Z + 0.5).astype(np.int16)
 
     # 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)
+    M = np.diag(imageS_mat.diagonal() / Z)    
+
+    image = np.zeros((D[2],D[1],D[0]),np.uint8)
     
     mode  = 2
     scale = 0
@@ -121,11 +122,9 @@
     # use the 6 dim parm_vector (3 angles, 3 translations) to remap
     #
     M_inverse = get_inverse_mappings(parm_vector)
-    (layers, rows, cols) = image.shape
 
     # allocate the zero image
-    #remaped_image = np.zeros(image.size, dtype=np.uint8).reshape(layers, rows, cols)
-    remaped_image = np.empty(image.shape, dtype=np.uint8)
+    remaped_image = np.zeros(image.shape, dtype=np.uint8)
 
     step = np.array([1, 1, 1], dtype=np.int32)
 
@@ -174,18 +173,18 @@
     M_inverse = build_rotate_matrix(-parm_vector)
     return M_inverse
 
-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'):
+def register(image1, image1_mat, image2, image2_mat, multires=[4, 2], histo_fwhm=3, 
+             ftype=1, lite=0, smhist=0, method='nmi', opt_method='hybrid',
+	     optimize_function=None):
 
     """
-    parm_vector = coregister(image1, image1_mat, image2, image2_mat, multires=[4, 2], histo_fwhm=3,
+    parm_vector = register(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 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 
-    package. The optimal parameter is returned.
+    package. The optimal rigid body parameter is returned.
 
     Parameters 
     ----------
@@ -216,7 +215,7 @@
         information; mi is mutual information; ecc is entropy cross
         correlation; ncc is normalized cross correlation. mse is mean
 	squared error.
-    opt_method: {'powell', 'hybrid'}, optional
+    opt_method: {'powell', 'cg', '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
         2 are powell, in hybrid pass 2 is conjugate gradient.
@@ -235,65 +234,124 @@
     >>> import _registration as reg
 
     >>> image1, image2, fwhm, improc = reg.demo_build_dual_volumes()
-    >>> parm_vector = coregister(image1, image2, fwhm, improc)
+    >>> parm_vector = register(image1, image2, fwhm, improc)
 
     """
 
-    start = time.time()
-    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)
+    # do the parameter validity checking. this is specific to this 3D registration.
+    # make sure the image is 3D and the mats are 4x4 with nonzero diagonal
+
+    if image1.ndim != 3:
+        raise ValueError, "Image 1 is not 3 dimensional"
+
+    if image2.ndim != 3:
+        raise ValueError, "Image 2 is not 3 dimensional"
+
+    if image1.dtype != np.uint8:
+        raise ValueError, "Image 1 is not 8 bit (required for joint histogram)"
+
+    if image2.dtype != np.uint8:
+        raise ValueError, "Image 2 is not 8 bit (required for joint histogram)"
+
+    if image1_mat.shape != (4,4):
+        raise ValueError, "Image1 MAT is not 4x4"
+
+    if image2_mat.shape != (4,4):
+        raise ValueError, "Image2 MAT is not 4x4"
+
+    if (np.diag(image1_mat)).prod() == 0:
+        raise ValueError, "Image1 MAT has a 0 on the diagonal"
+
+    if (np.diag(image2_mat)).prod() == 0:
+        raise ValueError, "Image2 MAT has a 0 on the diagonal"
+
+    if opt_method=='hybrid' and np.size(multires) != 2:
+        raise ValueError, "hybrid method must be 2 pass registration"
+
+    if ftype != 0 and ftype != 1: 
+        raise ValueError, "choose filter type 0 or 1 only"
+
+    if lite != 0 and lite != 1: 
+        raise ValueError, "choose histogram generation type 0 or 1 only"
+
+    if smhist != 0 and smhist != 1: 
+        raise ValueError, "choose histogram smoothing type 0 or 1 only"
+
+    if method != 'nmi' and method != 'mi'  and method != 'ncc'\
+                       and method != 'ecc' and method != 'mse':
+        raise ValueError, "choose cost method nmi, mi, ecc, mse, ncc"
+
+    if opt_method != 'powell' and opt_method != 'cg'  and opt_method != 'hybrid':
+        raise ValueError, "only optimize methods powell, cg or hybrid are supported"
+
+    # default is to use the cost_function I provided.
+    # this shows you can override this but the parameters will have to
+    # be changed for the new cost function if it is different
+
+    if optimize_function == None:
+        optimize_function = cost_function
+
+    parm_vector = multires_registration(optimize_function, image1, image1_mat, image2, image2_mat,
+		                        multires, histo_fwhm, lite, smhist, method, opt_method)
+
     return parm_vector
 
-def multires_registration(image1, image1_mat, image2, image2_mat, multires, histo_fwhm, 
-		          lite, smhist, method, opt_method):
+def multires_registration(optimize_function, 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 coregister() which optionally does 3D image filtering and 
-    provies timing for registration.
+    to be called by register() which does parameter validation 
 
     Parameters 
     ----------
-
     image1 : {nd_array} 
         image1 is the source image to be remapped during the registration. 
+    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.
-    lite : {integer}
+    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.
+    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: {integer}
+    smhist: {0, 1}, optional
         flag for joint histogram low pass filtering. 0 for no filter,
         1 for do filter.
-    method: {'nmi', 'mi', 'ncc', 'ecc', 'mse'}
+    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. mse is mean
-	square error.
-    opt_method: {'powell', 'hybrid'}
+	squared error.
+    opt_method: {'powell', 'cg', '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
         2 are powell, in hybrid pass 2 is conjugate gradient.
 
+
     Returns 
     -------
-    x : {nd_array}
+    parm_vector : {nd_array}
         this is the optimal alignment (6-dim) array with 3 angles and
         3 translations.
 
     Examples
     --------
 
-    (calling this from coregister which optionally filters image2)
+    (calling this from register which optionally filters image2)
     >>> import numpy as NP
     >>> import _registration as reg
     >>> image1, mat1, image2, mat2 = reg.demo_build_dual_volumes()
-    >>> parm_vector = coregister(image1, image2, imdata)
+    >>> parm_vector = register(image1, image2, imdata)
 
     """
     ret_histo=0
@@ -308,8 +366,10 @@
     for i in loop:
 	# this is the volume subsample
 	step[:] = multires[i]
-        optfunc_args = (image1, image1_mat, image2, image2_mat, step, fwhm, lite,
-			smhist, method, ret_histo)
+	# optfunc_args is specific to the cost_function in this file
+	# this will need to change if you use another optimize_function.
+        optfunc_args = (image1, image1_mat, image2, image2_mat, step, histo_fwhm,
+			lite, smhist, method, ret_histo)
         p_args = (optfunc_args,)
         if opt_method=='powell':
             print 'POWELL multi-res registration step size ', step
@@ -324,8 +384,8 @@
                 print 'Hybrid POWELL multi-res registration step size ', step
                 print 'vector ', x
                 lite = 0
-                optfunc_args = (image1, image1_mat, image2, image2_mat, step, 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_powell(optimize_function, x, args=p_args, callback=callback_powell) 
             elif i==1:
@@ -417,7 +477,7 @@
 
     return kernel
 
-def filter_image_3D(imageRaw, fwhm, ftype=2):
+def filter_image_3D(imageRaw, fwhm, ftype=2, give_2D=0):
     """
     image_F_xyz = filter_image_3D(imageRaw, fwhm, ftype=2):
     does 3D separable digital filtering using scipy.ndimage.correlate1d
@@ -427,9 +487,9 @@
     imageRaw : {nd_array}
         the unfiltered 3D volume image
     fwhm : {int}
-        used for kernel width
+        used for kernel width. this is 3 elements (one for each dimension)
     ktype: {1, 2}, optional
-        kernel type. 1 is Gauss convoled with spline, 2 is Gauss
+        kernel type. 1 is Gauss convoled with spline (SPM), 2 is Gauss
 
     Returns 
     -------
@@ -442,32 +502,37 @@
     >>> import _registration as reg
     >>> image1, image2, imdata = reg.demo_build_dual_volumes()
     >>> ftype = 1
-    >>> image_Filter_xyz = filter_image_3D(image1['data'], image1['fwhm'], ftype)
+    >>> image_Filter_xyz = filter_image_3D(image, fwhm, ftype)
     >>> image1['data'] = image_Filter_xyz
     """
 
-    p = np.ceil(2*fwhm[0]).astype(int)
-    x = np.array(range(-p, p+1))
+    p = np.ceil(2*fwhm).astype(int)
+    x = np.array(range(-p[0], p[0]+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))
+    x = np.array(range(-p[1], p[1]+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))
+    x = np.array(range(-p[2], p[2]+1))
     kernel_z = smooth_kernel(fwhm[2], x, ktype=ftype)
 
     output=None
-    # 3D filter in 3 1D separable stages
+    # 3D filter in 3 1D separable stages. keep the image
+    # names at each stage separate in case you need them
+    # for example may need an image that is 2D slice filtered only
     axis = 0
     image_F_x   = correlate1d(imageRaw,   kernel_x, axis, output)
     axis = 1
     image_F_xy  = correlate1d(image_F_x,  kernel_y, axis, output)
     axis = 2
     image_F_xyz = correlate1d(image_F_xy, kernel_z, axis, output)
-    return image_F_xyz  
 
+    if give_2D==0:
+        return image_F_xyz  
+    else:
+        return image_F_xyz, image_F_xy
+
+
 def build_fwhm(M, S):
     """
     fwhm = build_fwhm(M, S)
@@ -512,10 +577,10 @@
     # return the 3D Gaussian kernel width (xyz)
     return fwhm 
 
-def optimize_function(x, optfunc_args):
+def cost_function(x, optfunc_args):
     """
-    cost = optimize_function(x, optfunc_args)    --- OR ---
-    cost, joint_histogram = optimize_function(x, optfunc_args)   
+    cost = cost_function(x, optfunc_args)    --- OR ---
+    cost, joint_histogram = cost_function(x, optfunc_args)   
 
     computes the alignment between 2 volumes using cross correlation or mutual
     information metrics. In both the 8 bit joint histogram of the 2 images is
@@ -590,7 +655,7 @@
     >>> ret_histo = 1
     >>> optfunc_args = (image1, image2, imdata['step'], imdata['fwhm'], lite, smhist, method, ret_histo)
     >>> x = np.zeros(6, dtype=np.float64)
-    >>> return cost, joint_histogram = reg.optimize_function(x, optfunc_args)
+    >>> return cost, joint_histogram = reg.cost_function(x, optfunc_args)
 
 
     """
@@ -625,8 +690,7 @@
 
         # allocate the zero image
         #(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)
+        remap_image_F = np.zeros(image_F.shape, dtype=np.uint8)
         # trilinear interpolation mapping.
         reg.register_linear_resample(image_F, remap_image_F, composite, sample_vector)
         cost = (np.square(image_G-remap_image_F)).mean()
@@ -792,7 +856,7 @@
     return rot_matrix
 
 
-def build_test_volume(imagedesc, S=[15.0, 25.0, 10.0]):
+def build_test_volume(imagedesc, S=[1500.0, 2500.0, 1000.0]):
 
     """
     build a 3D Gaussian volume. user passes in image dims in imagedesc
@@ -831,21 +895,17 @@
     aa       = (np.square(a))/sigma[0]
     bb       = (np.square(b))/sigma[1]
     cc       = (np.square(c))/sigma[2]
-    volume3D = np.exp(-(aa + bb + cc))
+    volume3D = (255.0*np.exp(-(aa + bb + cc))).astype(np.uint8)
 
     return volume3D
 
 
 
-def load_volume(imagedesc, imagename=None, threshold=0.999, debug=0):
+def load_volume(imagedesc, imagename=None):
 
     """
-    image = load_volume(imagedesc, imagename=None, threshold=0.999, debug=0)  --- OR ---
-    image, h, ih, index = load_volume(imagedesc, imagename=None, threshold=0.999, debug=0)
 
-    gets an image descriptor and optional filename and returns a scaled 8 bit volume. The
-    scaling is designed to make full use of the 8 bits (ignoring high amplitude outliers).
-    The current method uses numpy fromfile and will be replaced by neuroimage nifti load.
+    returns an nd_array from the filename or blank image. used for testing. 
 
     Parameters 
     ----------
@@ -856,47 +916,21 @@
         name of image file. No name creates a blank image that is used for creating
         a rotated test image or image rescaling.
 
-    threshold : {float} : optional
-        this is the threshold for upper cutoff in the 8 bit scaling. The volume histogram
-        and integrated histogram is computed and the upper amplitude cutoff is where the 
-        integrated histogram crosses the value set in the threshold. setting threshold to
-        1.0 means the scaling is done over the min to max amplitude range.
-
-    debug : {0, 1} : optional
-        when debug=1 the method returns the volume histogram, integrated histogram and the 
-        amplitude index where the provided threshold occured.
-
     Returns 
     -------
     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)
+    M : {nd_array}
+        the voxel-to-physical affine matrix (mat)
 
-    image : {nd_array}
-        the volume data assoicated with the filename or a blank volume of the same
-        dimensions as specified in imagedesc.
-
-    h : {nd_array}
-        the volume 1D amplitude histogram
-
-    ih : {nd_array}
-        the volume 1D amplitude integrated histogram
-
-    index : {int}
-        the amplitude (histogram index) where the integrated histogram
-        crosses the 'threshold' provided.
-
     Examples
     --------
 
-    >>> import numpy as NP
     >>> import _registration as reg
     >>> anat_desc = reg.load_anatMRI_desc()
-    >>> image_anat, h, ih, index = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img', debug=1)
-    >>> index
-    210
+    >>> image, M = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img')
 
 
     """
@@ -906,8 +940,6 @@
     if imagename == None:
         # imagename of none means to create a blank image
         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:
         image = np.fromfile(imagename,
                    dtype=np.uint16).reshape(imagedesc['layers'], imagedesc['rows'], imagedesc['cols']);
@@ -922,32 +954,133 @@
     if imagename == None:
         # no voxels to scale to 8 bits
         image = image.astype(np.uint8)
-        return image, M
 
-    # 8 bit scale with threshold clip of the volume integrated histogram
+    return image, M
+
+
+
+def scale_image(image, max_amp=255, image_type=np.uint8, threshold=0.999, fetch_ih=0):
+
+    """
+    scale and threshold clip the volume using the integrated histogram
+    to set the high threshold
+
+    Parameters 
+    ----------
+    image : {nd_array}
+        raw unscaled volume
+
+    max_amp : int (default 255)
+        the maximum value of the scaled image
+
+    image_type : nd_array dtype (default uint8)
+        the type of the volume to return.
+
+    threshold : float (default 0.999)
+        the value of the normalized integrated histogram
+	that when reached sets the high threshold index
+
+    Returns 
+    -------
+    image : {nd_array}
+        the scaled volume
+    ih : {nd_array}
+        the integrated histogram. can be used for image display 
+	purpose (histogram equalization)
+
+    """
+
     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:
-        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  = reg.register_image_threshold(image, h, ih, threshold)
-    scale  = 255.0 / (index-min)
-    # generate the scaled 8 bit image
-    image = (scale*(image.astype(np.float)-min))
-    image[image>255] = 255 
-    image = image.astype(np.uint8)
-    if debug == 1:
-        return image, M, h, ih, index
+    if max == 0 and min == 0:
+        raise ValueError, "Zero image. cannot be scaled"
+
+    # need range of pixels for the number of bins
+    h, edges = np.histogram(image, bins=(max-min))
+    ih = (np.cumsum(h)).astype(np.float64)
+    # normalize the integrated histogram
+    ih = ih / ih.max()
+    indices = np.where(ih >= threshold)
+    # wind up getting all the indices where the ih >= threshold
+    # and only need the first index. tuple has one nd_array and
+    # get the 0 element from it ([0][0])
+    index   = indices[0][0]
+    scale   = float(max_amp) / (index-min)
+    image   = (scale*(image.astype(np.float)-min))
+    image[image>max_amp] = max_amp
+    # down type. usually will go from float to 8 bit (needed for the 8 bit joint histogram)
+    image = image.astype(image_type)
+
+    if fetch_ih == 1:
+        return image, ih
     else:
-        return image, M
+        return image
 
 
+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):
+                    
+    """
+    test the cost function and (optional) view the joint histogram. can be used
+    during intra-modal registration to measure the current alignment (return
+    the cross correlation). would measure before and after registration
 
+
+
+    """
+
+    # do the parameter validity checking. this is specific to this 3D registration.
+    # make sure the image is 3D and the mats are 4x4 with nonzero diagonal
+
+    if image1.ndim != 3:
+        raise ValueError, "Image 1 is not 3 dimensional"
+
+    if image2.ndim != 3:
+        raise ValueError, "Image 2 is not 3 dimensional"
+
+    if image1.dtype != np.uint8:
+        raise ValueError, "Image 1 is not 8 bit (required for joint histogram)"
+
+    if image2.dtype != np.uint8:
+        raise ValueError, "Image 2 is not 8 bit (required for joint histogram)"
+
+    if image1_mat.shape != (4,4):
+        raise ValueError, "Image1 MAT is not 4x4"
+
+    if image2_mat.shape != (4,4):
+        raise ValueError, "Image2 MAT is not 4x4"
+
+    if (np.diag(image1_mat)).prod() == 0:
+        raise ValueError, "Image1 MAT has a 0 on the diagonal"
+
+    if (np.diag(image2_mat)).prod() == 0:
+        raise ValueError, "Image2 MAT has a 0 on the diagonal"
+
+    if method != 'nmi' and method != 'mi'  and method != 'ncc'\
+                       and method != 'ecc' and method != 'mse':
+        raise ValueError, "choose cost method nmi, mi, ecc, mse, ncc"
+
+    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)
+			
+    if ret_histo:
+        cost, joint_histogram = cost_function(P, optfunc_args)
+        return cost, joint_histogram 
+    else:
+        cost = cost_function(P, optfunc_args)
+        return cost
+
+
+
 #
 #  ---- demo/debug routines  ---- 
 #
@@ -981,50 +1114,20 @@
     return files_fMRI
 
 
-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 (optional) view the joint histogram
-    # default of use of ncc for testing the cross-correlation as a metric
-    # of alignment
-    #
 
-    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(P, optfunc_args)
-        return cost, joint_histogram 
-    else:
-        cost = optimize_function(P, optfunc_args)
-        return cost
-
 def build_scale_volume(image, mat, scale):
     #
     # rescale the 'mat' (voxel to physical mapping matrix) 
     #
+    M = mat * scale
     (layers, rows, cols) = image.shape
-    M = mat * scale
     # dimensions 
     D = np.zeros(3, dtype=np.int32);
     Z = np.zeros(3, dtype=np.float64);
     D[0] = rows/scale
     D[1] = cols/scale
     D[2] = layers/scale
-    #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)
+    image2 = np.zeros([D[2], D[0], D[1]], dtype=np.uint8)
     mode = 1;
     reg.register_volume_resample(image, image2, Z, scale, mode)
     return image2, M
@@ -1036,7 +1139,7 @@
     builds a volume and a scaled-rotated version for coreg testing 
 
     image1, mat1, image2, mat2 = reg.demo_build_dual_volumes()
-    x = reg.coregister(image1, mat1, image2, mat2, method='ncc', lite=1) 
+    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)
 
@@ -1050,6 +1153,7 @@
     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
@@ -1064,12 +1168,75 @@
     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
@@ -1142,8 +1309,12 @@
 
     # 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)
@@ -1174,7 +1345,7 @@
         # 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 = coregister(imageF, fmri_mat, imageG, fmri_mat, lite=histo_method, method='ncc',
+        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,
@@ -1200,7 +1371,7 @@
     if 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,
+    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 '



More information about the Scipy-svn mailing list