[Scipy-svn] r4638 - in trunk/scipy/ndimage: . src tests

scipy-svn@scip... scipy-svn@scip...
Mon Aug 11 18:09:06 CDT 2008


Author: jarrod.millman
Date: 2008-08-11 18:08:56 -0500 (Mon, 11 Aug 2008)
New Revision: 4638

Removed:
   trunk/scipy/ndimage/_registration.py
   trunk/scipy/ndimage/_segmenter.py
   trunk/scipy/ndimage/src/register/
   trunk/scipy/ndimage/src/segment/
   trunk/scipy/ndimage/tests/test_registration.py
   trunk/scipy/ndimage/tests/test_segment.py
Modified:
   trunk/scipy/ndimage/setup.py
Log:
removing segmentation and registration code since it isn't ready for release.


Deleted: trunk/scipy/ndimage/_registration.py
===================================================================
--- trunk/scipy/ndimage/_registration.py	2008-08-11 21:39:17 UTC (rev 4637)
+++ trunk/scipy/ndimage/_registration.py	2008-08-11 23:08:56 UTC (rev 4638)
@@ -1,1044 +0,0 @@
-#
-# written by Tom Waite
-# rigid body 3D registration
-#
-
-
-import math
-import numpy as np
-from scipy.special import erf
-from scipy.ndimage import correlate1d
-from scipy.optimize import fmin_powell, fmin_cg
-
-import scipy.ndimage._register as reg
-
-import time
-import glob
-
-# Issue warning regarding heavy development status of this module
-import warnings
-_msg = "The registration code is under heavy development and therefore the \
-public API will change in the future.  The NIPY group is actively working on \
-this code, and has every intention of generalizing this for the Scipy \
-community.  Use this module minimally, if at all, until it this warning is \
-removed."
-warnings.warn(_msg, UserWarning)
-
-# TODO:  Add docstrings for public functions in extension code.
-# Add docstrings to extension code.
-#from numpy.lib import add_newdoc
-#add_newdoc('scipy.ndimage._register', 'register_histogram',
-#    """A joint histogram used for registration module.
-#    """)
-
-
-#
-#  ---- co-registration and IO  ---- 
-#
-
-def resize_image(imageS, imageS_mat, imageR_mat):
-    """
-    zoom_image = resize_image(imageS, imageS_mat, imageR_mat)
-
-    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 
-    ----------
-
-    imageS: {ndarray} 
-        imageS is the source image to be resized.
-
-    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 : {ndarray}
-
-    Examples
-    --------
-
-    >>> import _registration as reg
-    >>> measures, image_anat, image_anat_mat, image_fmri_mat, fmri_series = reg.demo_MRI_coregistration()
-
-    >>> resampled_fmri = reg.resize_image(fmri_series[10], image_fmri_mat, image_anat_mat)
-
-    image 10 in the fmri_series is resampled from image_fmri_mat to image_anat coordinates
-
-    """
-
-    # get the zoom
-    Z = imageS_mat.diagonal() / imageR_mat.diagonal()
-
-    # 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)
-
-    # for the test data, set the xyz voxel sizes for fMRI volume. M is a 4x4 matrix.
-    M = np.diag(imageS_mat.diagonal() / Z)    
-
-    image = np.zeros((D[2],D[1],D[0]),np.uint8)
-    
-    mode  = 2
-    scale = 0
-    reg.register_volume_resample(imageS, image, Z, scale, mode)
-
-    return image, M
-
-def remap_image(image, M_inverse, resample='linear'):
-    """
-    remaped_image = remap_image(image, parm_vector, resample='linear')
-
-    rotates and translates image using the 3 angles and 3 translations in the 6-dim
-    parm_vector. The mapping is stored in the 4x4 M_inverse matrix from the get_inverse_mapping
-    method.
-
-    Parameters 
-    ----------
-    image : {ndarray} 
-        image is the source image to be remapped. 
-
-    M_inverse : {ndarray}
-        M_inverse is the 4x4 inverse affine matrix 
-
-    resample : {'linear', 'cubic'}, optional
-
-
-    Returns 
-    -------
-    remaped_image : {ndarray}
-
-    Examples
-    --------
-        image = fmri_series[i]
-        x[0:6] = measures[i]['align_rotate'][0:6]
-        M = get_inverse_mappings(x)
-        # overwrite the fMRI volume with the aligned volume
-        fmri_series[i] = remap_image(image, M, resample='cubic')
-
-    """
-
-    # allocate the zero image
-    remaped_image = np.zeros(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, remaped_image, M_inverse, step)
-    elif resample == 'cubic':
-        # tricubic convolve interpolation mapping. 
-        reg.register_cubic_resample(image, remaped_image, M_inverse, step)
-
-    return remaped_image
-
-
-def get_inverse_mappings(parm_vector):
-    """
-    M_inverse = get_inverse_mappings(parm_vector)
-
-    takes the 6-dimensional rotation/translation vector and builds the inverse
-    4x4 mapping matrix M_inverse that will map imageG to imageF orientation
-
-    Parameters 
-    ----------
-    parm_vector : {nd_array} 
-
-    Returns 
-    -------
-    M_inverse : {nd_array}
-
-    Examples
-    --------
-
-    >>> import numpy as NP
-    >>> import _registration as reg
-    >>> array = np.zeros(6, dtype=float)
-    >>> M = reg.get_inverse_mappings(array)
-    >>> M 
-
-    array([
-    [ 1.,  0.,  0.,  0.],
-    [ 0.,  1.,  0.,  0.],
-    [ 0.,  0.,  1.,  0.],
-    [ 0.,  0.,  0.,  1.]])
-
-    """
-    # get the inverse mapping to rotate the G matrix to F space following registration
-    # -parm_vector is the inverse angles and translations
-    M_inverse = build_rotate_matrix(-parm_vector)
-    return M_inverse
-
-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 = register(image1, image1_mat, image2, image2_mat, multires=[4, 2], histo_fwhm=3,
-                             ftype=1, lite=0, smhist=0, method='nmi', opt_method='powell'):
-
-    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 rigid body parameter is returned.
-
-    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. 
-    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: {0, 1}, optional
-        flag for joint histogram low pass filtering. 0 for no filter,
-        1 for do filter.
-    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
-        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 
-    -------
-    parm_vector : {nd_array}
-        this is the optimal alignment (6-dim) array with 3 angles and
-        3 translations.
-
-    Examples
-    --------
-
-    >>> import numpy as NP
-    >>> import _registration as reg
-
-    >>> image1, image2, fwhm, improc = reg.demo_build_dual_volumes()
-    >>> parm_vector = register(image1, image2, fwhm, improc)
-
-    """
-
-    # 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(optimize_function, image1, image1_mat, image2, image2_mat,
-                          multires, histo_fwhm, lite, smhist, method, opt_method):
-
-    """
-
-    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. 
-    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: {0, 1}, optional
-        flag for joint histogram low pass filtering. 0 for no filter,
-        1 for do filter.
-    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
-        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 
-    -------
-    parm_vector : {nd_array}
-        this is the optimal alignment (6-dim) array with 3 angles and
-        3 translations.
-
-    Examples
-    --------
-
-    (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 = register(image1, image2, imdata)
-
-    """
-    ret_histo=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(np.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:
-        # this is the volume subsample
-        step[:] = multires[i]
-        # 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
-            print 'vector ', x
-            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
-            x = fmin_cg(optimize_function, x, args=p_args, callback=callback_cg) 
-        elif opt_method=='hybrid':
-            if i==0:
-                print 'Hybrid POWELL multi-res registration step size ', step
-                print 'vector ', x
-                lite = 0
-                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:
-                print 'Hybrid CG multi-res registration step size ', step
-                print 'vector ', x
-                lite = 1
-                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) 
-
-    return x
-
-
-def callback_powell(x):
-    """
-    called by optimize.powell only. prints current parameter vector.
-    """
-    print 'Parameter Vector from Powell: - '
-    print x
-    return
-
-def callback_cg(x):
-    """
-    called by optimize.cg only. prints current parameter vector.
-    """
-    print 'Parameter Vector from Conjugate Gradient: - '
-    print x
-    return
-
-def smooth_kernel(fwhm, x, pixel_scale=8.0, ktype=1):
-    """
-    kernel = smooth_kernel(fwhm, x, ktype=1)
-
-    smooth_kernel creates filter kernel based on image sampling parameters.
-    provide domain of kernel and sampling parameters. 
-
-    Parameters 
-    ----------
-    fwhm : {int}
-        used for kernel width
-    x : {nd_array}
-        domain of kernel
-    ktype: {1, 2}, optional
-        kernel type. 1 is Gauss convoled with spline, 2 is Gauss
-
-
-    Returns 
-    -------
-    kernel : {nd_array}
-
-    Examples
-    --------
-
-    >>> import numpy as NP
-    >>> import _registration as reg
-    >>> fwhm = 3
-    >>> ftype = 2
-    >>> p = np.ceil(2*fwhm).astype(int)
-    >>> x = np.array(range(-p, p+1))
-    >>> kernel = reg.smooth_kernel(fwhm, x, ktype=ftype)
-    >>> kernel
-
-    array([
-     4.77853772e-06,   1.41575516e-04,   2.26516955e-03,
-     1.95718875e-02,   9.13238336e-02,   2.30120330e-01,
-     3.13144850e-01,   2.30120330e-01,   9.13238336e-02,
-     1.95718875e-02,   2.26516955e-03,   1.41575516e-04,
-     4.77853772e-06])
-
-    """
-    eps = 0.00001
-    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)
-        w2 = -0.5 / s
-        w3 = math.sqrt((s*math.pi) /2.0)
-        kernel = 0.5*(erf(w1*(x+1))*(x+1) + erf(w1*(x-1))*(x-1)
-                      - 2.0*erf(w1*(x))*(x) + w3*(np.exp(w2*np.square(x+1))) 
-                      + np.exp(w2*(np.square(x-1)))
-                      - 2.0*np.exp(w2*np.square(x)))
-        kernel[kernel<0] = 0
-        kernel = kernel / kernel.sum()  
-    else:
-        # Gauss kernel 
-        kernel = (1.0/math.sqrt(2.0*math.pi*s)) * np.exp(-np.square(x)/(2.0*s)) 
-        kernel = kernel / kernel.sum()  
-
-    return kernel
-
-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
-
-    Parameters 
-    ----------
-    imageRaw : {nd_array}
-        the unfiltered 3D volume image
-    fwhm : {int}
-        used for kernel width. this is 3 elements (one for each dimension)
-    ktype: {1, 2}, optional
-        kernel type. 1 is Gauss convoled with spline (SPM), 2 is Gauss
-
-    Returns 
-    -------
-    image_F_xyz : {nd_array}
-        3D filtered volume image
-
-    Examples
-    --------
-
-    >>> import _registration as reg
-    >>> image1, image2, imdata = reg.demo_build_dual_volumes()
-    >>> ftype = 1
-    >>> image_Filter_xyz = filter_image_3D(image, fwhm, ftype)
-    >>> image1['data'] = image_Filter_xyz
-    """
-
-    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)
-
-    x = np.array(range(-p[1], p[1]+1))
-    kernel_y = smooth_kernel(fwhm[1], x, ktype=ftype)
-
-    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. 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)
-
-    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)
-
-    builds the low pass filter kernel sigma value from the image pixel sampling
-
-    Parameters 
-    ----------
-    M : {nd_array}
-        input 4x4 voxel to physical map matrix (called 'MAT')
-
-    S : {nd_array}
-        1x3 sample increment array. should be = (1, 1, 1)
-
-    Returns 
-    -------
-    fwhm : {nd_array}
-        the 3D Gaussian kernel width
-
-    Examples
-    --------
-
-    >>> import numpy as NP
-    >>> import _registration as reg
-    >>> anat_desc = reg.load_anatMRI_desc()
-    >>> image1 = reg.load_volume(anat_desc, imagename='ANAT1_V0001.img')
-    >>> 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 in the first row
-    vxg = np.sqrt(view_3x3.sum(axis=1))
-    # 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
-    x[x<0] = 0
-    fwhm = np.sqrt(x) / vxg
-    # pathology when stepsize = 1 for MAT equal to the identity matrix
-    fwhm[fwhm==0] = 1
-    # return the 3D Gaussian kernel width (xyz)
-    return fwhm 
-
-def cost_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
-    computed. The 8 bit scaling is done using an integrated histogram method and
-    is called prior to this.
-
-    Parameters 
-    ----------
-    x : {nd_array}
-        this is the current (6-dim) array with 3 angles and 3 translations.
-
-    optfunc_args : {tuple}
-        this is a tuple of 8 elements that is formed in the scipy.optimize powell
-        and cg (conjugate gradient) functions. this is why the elements need to be
-        a tuple. The elements of optfunc_args are: 
-
-        image_F       : {dictionary} 
-            image_F is the source image to be remapped during the registration. 
-            it is a dictionary with the data as an ndarray in the ['data'] component.
-        image_G       : {dictionary} 
-            image_G is the reference image that image_F gets mapped to. 
-        sample_vector : {nd_array} 
-            sample in x,y,x. should be (1,1,1)
-        fwhm          : {nd_array} 
-            Gaussian sigma
-        do_lite       : {0, 1} 
-            lite of 1 is to jitter both images during resampling. 
-            0 is to not jitter. jittering is for non-aliased volumes.
-        smooth        : {0, 1} 
-            flag for joint histogram low pass filtering. 0 for no filter,
-            1 for do filter.
-        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. mse is mean
-            square error. with mse there is no joint histogram.
-        ret_histo     : {0, 1} 
-            if 0 return is: cost 
-            if 0 return is: cost, joint_histogram  
-
-    Returns 
-    -------
-        cost : {float}
-            the negative of one of the mutual information metrics
-            or negative cross correlation. use negative as the optimization
-            is minimization.
-
-        --- OR --- (if ret_histo = 1)
-
-        cost : {float}
-            the negative of one of the mutual information metrics
-            or negative cross correlation. use negative as the optimization
-            is minimization.
-
-        joint_histogram : {nd_array}
-            the 2D (256x256) joint histogram of the two volumes
-
-
-    Examples
-    --------
-
-    >>> import numpy as NP
-    >>> import _registration as reg
-    >>> 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')
-    >>> image1['fwhm'] = reg.build_fwhm(image1['mat'], imdata['step'])
-    >>> image2['fwhm'] = reg.build_fwhm(image2['mat'], imdata['step'])
-    >>> method = 'ncc'
-    >>> lite = 1
-    >>> smhist = 0
-    >>> 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.cost_function(x, optfunc_args)
-
-
-    """
-
-    image_F       = optfunc_args[0]
-    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 (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)
-    composite = np.dot(composite, rot_matrix)
-
-    if method == 'mse':
-        #
-        # mean squard error method
-        #
-
-        # allocate the zero image
-        #(layers, rows, cols) = image_F.shape
-        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()
-        # cost is min when G and F are aligned so keep cost positive
-
-        return cost
-
-    else:
-        #
-        # histogram-based methods (nmi, ncc, mi, ecc)
-        #
-
-        # allocate memory for 2D histogram
-        joint_histogram = np.zeros([256, 256], dtype=np.float64)
-
-        if do_lite: 
-            reg.register_histogram_lite(image_F, image_G, composite, sample_vector, joint_histogram)
-        else:
-            reg.register_histogram(image_F, image_G, composite, sample_vector, joint_histogram)
-
-        # smooth the histogram
-        if smooth: 
-            p = np.ceil(2*fwhm).astype(int)
-            x = np.array(range(-p, p+1))
-            hkernel = smooth_kernel(fwhm, x)
-            output=None
-            # 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
-            joint_histogram = correlate1d(joint_histogram, hkernel, axis, output)
-            axis = 1
-            joint_histogram = correlate1d(joint_histogram, hkernel, axis, output)
-
-        joint_histogram += epsilon # prevent log(0) 
-        # normalize the joint histogram
-        joint_histogram /= joint_histogram.sum() 
-        # get the marginals
-        marginal_col = joint_histogram.sum(axis=0)
-        marginal_row = joint_histogram.sum(axis=1)
-
-        if method == 'mi':
-            # mutual information
-            marginal_outer = np.outer(marginal_col, marginal_row)
-            H = joint_histogram * np.log(joint_histogram / marginal_outer)  
-            mutual_information = H.sum()
-            cost = -mutual_information
-
-        elif method == 'ecc':
-            # entropy correlation coefficient 
-            marginal_outer = np.outer(marginal_col, marginal_row)
-            H = joint_histogram * np.log(joint_histogram / marginal_outer)  
-            mutual_information = H.sum()
-            row_entropy = marginal_row * np.log(marginal_row)
-            col_entropy = marginal_col * np.log(marginal_col)
-            ecc  = -2.0*mutual_information/(row_entropy.sum() + col_entropy.sum())
-            cost = -ecc
-
-        elif method == 'nmi':
-            # normalized mutual information
-            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())
-            cost = -nmi
-
-        elif method == 'ncc':
-            # cross correlation from the joint histogram 
-            r, c = joint_histogram.shape
-            i = np.array(range(1,c+1))
-            j = np.array(range(1,r+1))
-            m1 = (marginal_row * i).sum()
-            m2 = (marginal_col * j).sum()
-            sig1 = np.sqrt((marginal_row*(np.square(i-m1))).sum())
-            sig2 = np.sqrt((marginal_col*(np.square(j-m2))).sum())
-            [a, b] = np.mgrid[1:c+1, 1:r+1]
-            a = a - m1
-            b = b - m2
-            # element multiplies in the joint histogram and grids
-            H = ((joint_histogram * a) * b).sum()
-            ncc  = H / (np.dot(sig1, sig2)) 
-            cost = -ncc
-
-        if ret_histo:
-            return cost, joint_histogram 
-        else:
-            return cost
-
-
-def build_rotate_matrix(img_data_parms):
-    """
-    rot_matrix = reg.build_rotate_matrix(img_data_parms)
-
-    takes the 6 element vector (3 angles, 3 translations) and build the 4x4 mapping matrix
-
-    Parameters 
-    ----------
-    img_data_parms : {nd_array}
-        this is the current (6-dim) array with 3 angles and 3 translations.
-
-    Returns 
-    -------
-    rot_matrix: {nd_array}
-        the 4x4 mapping matrix
-
-    Examples
-    --------
-
-    >>> import numpy as NP
-    >>> import _registration as reg
-    >>> x = np.zeros(6, dtype=np.float64)
-    >>> M = reg.build_rotate_matrix(x)
-    >>> M 
-    array([[ 1.,  0.,  0.,  0.],
-           [ 0.,  1.,  0.,  0.],
-           [ 0.,  0.,  1.,  0.],
-           [ 0.,  0.,  0.,  1.]])
-
-
-    """
-
-    R1 = np.zeros([4,4], dtype=np.float64);
-    R2 = np.zeros([4,4], dtype=np.float64);
-    R3 = np.zeros([4,4], dtype=np.float64);
-    T  = np.eye(4, dtype=np.float64);
-
-    alpha = math.radians(img_data_parms[0])
-    beta  = math.radians(img_data_parms[1])
-    gamma = math.radians(img_data_parms[2])
-
-    R1[0][0] = 1.0
-    R1[1][1] = math.cos(alpha)
-    R1[1][2] = math.sin(alpha)
-    R1[2][1] = -math.sin(alpha)
-    R1[2][2] = math.cos(alpha)
-    R1[3][3] = 1.0
-
-    R2[0][0] = math.cos(beta)
-    R2[0][2] = math.sin(beta)
-    R2[1][1] = 1.0
-    R2[2][0] = -math.sin(beta)
-    R2[2][2] = math.cos(beta)
-    R2[3][3] = 1.0
-
-    R3[0][0] = math.cos(gamma)
-    R3[0][1] = math.sin(gamma)
-    R3[1][0] = -math.sin(gamma)
-    R3[1][1] = math.cos(gamma)
-    R3[2][2] = 1.0
-    R3[3][3] = 1.0
-
-    T[0][0] = 1.0
-    T[1][1] = 1.0
-    T[2][2] = 1.0
-    T[3][3] = 1.0
-    T[0][3] = img_data_parms[3]
-    T[1][3] = img_data_parms[4]
-    T[2][3] = img_data_parms[5]
-
-    rot_matrix = np.dot(T, R1);
-    rot_matrix = np.dot(rot_matrix, R2);
-    rot_matrix = np.dot(rot_matrix, R3);
-
-    return rot_matrix
-
-
-def build_gauss_volume(imagedesc, S=[1500.0, 2500.0, 1000.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 = (255.0*np.exp(-(aa + bb + cc))).astype(np.uint8)
-
-    return volume3D
-
-
-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()
-    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
-
-
-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
-
-
-
-def build_scale_volume(image, mat, scale):
-    #
-    # rescale the 'mat' (voxel to physical mapping matrix) 
-    #
-    M = mat * scale
-    (layers, rows, cols) = image.shape
-    # 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[2], D[0], D[1]], dtype=np.uint8)
-    mode = 1;
-    reg.register_volume_resample(image, image2, Z, scale, mode)
-    return image2, M
-
-
-
-

Deleted: trunk/scipy/ndimage/_segmenter.py
===================================================================
--- trunk/scipy/ndimage/_segmenter.py	2008-08-11 21:39:17 UTC (rev 4637)
+++ trunk/scipy/ndimage/_segmenter.py	2008-08-11 23:08:56 UTC (rev 4638)
@@ -1,1718 +0,0 @@
-import math
-import numpy as NP
-import scipy.ndimage._segment as S
-
-_objstruct = NP.dtype([('Left', 'i'),
-                       ('Right', 'i'),
-                       ('Top', 'i'),
-                       ('Bottom', 'i'),
-                       ('Front', 'i'),
-                       ('Back', 'i'),
-                       ('Label', 'i'),
-                       ('Mass', 'i'),
-                       ('cX', 'f'),
-                       ('cY', 'f'),
-                       ('cZ', 'f'),
-                       ('voxelMean', 'f'),
-                       ('voxelVar', 'f'),
-                       ('COM', 'f', 6),
-                       ('TEM', 'f', 21)]
-                      )
-
-
-# Issue warning regarding heavy development status of this module
-import warnings
-_msg = "The segmentation code is under heavy development and therefore the \
-public API will change in the future.  The NIPY group is actively working on \
-this code, and has every intention of generalizing this for the Scipy \
-community.  Use this module minimally, if at all, until it this warning is \
-removed."
-warnings.warn(_msg, UserWarning)
-
-def canny_hysteresis(magnitude, canny_stats):
-    """
-    edge_image = canny_hysteresis(magnitude, canny_stats)
-
-    hystereis stage of Canny filter
-
-    Parameters 
-    ----------
-
-    magnitude : {nd_array}
-        the output from the canny_nonmax_supress() method
-
-    canny_stats : {dictionary}
-        contains the low and high thesholds determined from canny_nonmax_supress()
-
-    Returns 
-    ----------
-    edge_image : {nd_array}
-        the labeled edge image that can be displayed and used for later processing
-
-    """
-    [rows, cols] = magnitude.shape
-    edge_image = NP.zeros(rows*cols, dtype=NP.int16).reshape(rows, cols)
-    S.canny_hysteresis(magnitude, edge_image, canny_stats['low'], canny_stats['high'])
-
-    return edge_image
-
-def canny_nonmax_supress(horz_DGFilter, vert_DGFilter, img_means, thres=0.5, 
-                         mode=1, canny_l=0.5, canny_h=0.8):
-    """
-    magnitude, canny_stats = canny_nonmax_supress(horz_DGFilter, vert_DGFilter, img_means,
-                                          thres=0.5, mode=1, canny_l=0.5, canny_h=0.8)
-
-    non-max supression stage of Canny filter
-
-    Parameters 
-    ----------
-
-    horz_DGFilter : {nd_array}
-        the horizonal filtered image using the derivative of Gaussian kernel filter.
-        this is the output of the canny_filter method
-
-    vert_DGFilter : {nd_array}
-        the vertical filtered image using the derivative of Gaussian kernel filter.
-        this is the output of the canny_filter method
-
-    img_means : {dictionary}
-        mean X and Y values of edge signals determined from canny_filter
-
-    thres : {float}, optional
-        edge threshold applied to x and y filtered means
-
-    mode : {0, 1}, optional
-        threshold based on histogram mean(0) or mode(1)
-
-    canny_low : {float}, optional
-        low threshold applied to magnitude filtered 
-    
-    canny_high : {float}, optional
-        high threshold applied to magnitude filtered 
-
-    Returns 
-    ----------
-
-    magnitude : {nd_array}
-        magnitude of X and Y filtered for critical samples
-
-    canny_stats : {dictionary}
-        mean, low and high to be used for hysteresis
-
-    """
-    [rows, cols] = horz_DGFilter.shape
-    magnitude = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
-    aveMag, canny_low, canny_high = S.canny_nonmax_supress(horz_DGFilter, vert_DGFilter,
-                                        magnitude, mode, img_means['x-dg']*thres,
-                                        img_means['y-dg']*thres, canny_l, canny_h)
-
-    canny_stats = {'mean' : aveMag, 'low' : canny_low, 'high' : canny_high} 
-
-    return magnitude, canny_stats
-
-def canny_filter(slice, dg_kernel):
-    """
-    horz_DGFilter, vert_DGFilter, img_means = canny_filter(slice, dg_kernel)
-
-    Canny derivative of Gaussian filter 
-    returns the X and Y filterd image
-
-    Parameters 
-    ----------
-
-    slice : {nd_array}
-        2D image array
-
-    dg_kernel : {dictionary}
-        derivative of Gaussian kernel from build_d_gauss_kernel()
-
-    Returns 
-    ----------
-
-    horz_DGFilter : {nd_array}
-        X filtered image 
-
-    vert_DGFilter : {nd_array}
-        Y filtered image 
-
-    img_means : {dictionary}
-
-
-    """
-    slice = slice.astype(NP.float64)
-    [rows, cols] = slice.shape
-    horz_DGFilter = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
-    vert_DGFilter = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
-    aveX, aveY = S.canny_filter(slice, horz_DGFilter, vert_DGFilter,
-                   dg_kernel['coefficients'], dg_kernel['kernelSize'])
-
-    img_means = {'x-dg' : aveX, 'y-dg' : aveY} 
-
-    return horz_DGFilter, vert_DGFilter, img_means
-
-
-def binary_edge(label_image, ROI):
-    """
-    binary_edge_image = binary_edge(label_image, ROI)
-
-    takes the ROI dictionary with the blob bounding boxes and generates 
-    the binary edge for each blob. The ROI mask is used to isolate
-    the binary edges for later use (e.g. contour tracing).
-
-    Parameters 
-    ----------
-
-    label_image : {nd_array}
-        an image with labeled regions from get_blobs() method
-
-    ROI : {dictionary}
-        Region of Interest structure that has blob bounding boxes
-
-    Returns 
-    ----------
-
-    binary_edge_image : {nd_array}
-        edge image for each ROI combined into a single image
-
-    """
-
-    [rows, cols]      = label_image.shape
-    binary_edge_image = NP.zeros(rows*cols, dtype=NP.uint16).reshape(rows, cols)
-    number_regions    = ROI.size
-    indices           = range(0, number_regions)
-    for i in indices:
-        left   = ROI[i]['Left']-2
-        right  = ROI[i]['Right']+2
-        bottom = ROI[i]['Bottom']-2
-        top    = ROI[i]['Top']+2
-        Label  = ROI[i]['Label']
-        if left < 0: 
-            left = 0
-        if bottom < 0: 
-            bottom = 0
-        if right > cols-1: 
-            right = cols-1
-        if top > rows-1: 
-            top = rows-1
-
-        roi_rows = top-bottom
-        roi_cols = right-left
-        label_region  = NP.zeros(roi_rows*roi_cols, dtype=NP.uint16).reshape(roi_rows, roi_cols)
-        input = NP.zeros(roi_rows*roi_cols, dtype=NP.uint16).reshape(roi_rows, roi_cols)
-        # load the labeled region 
-        label_region[0:roi_rows, 0:roi_cols][label_image[bottom:top, left:right]==Label] = 1 
-        S.binary_edge(label_region, input)
-        input[0:roi_rows,0:roi_cols][input[0:roi_rows,0:roi_cols]==1] = Label
-        binary_edge_image[bottom:top,left:right] = binary_edge_image[bottom:top,left:right] + \
-                                                               input[0:roi_rows,0:roi_cols] 
-
-    return binary_edge_image
-
-
-def roi_co_occurence(label_image, raw_image, ROI, distance=2, orientation=90, verbose=0):
-    """
-    roi_co_occurence(label_image, raw_image, ROI, distance=2, orientation=90, verbose=0)
-
-    - OR -
-
-    texture_arrays = roi_co_occurence(label_image, raw_image, ROI, distance=2, verbose=1)
-
-    (N-S, E-W, NW-SE, NE-SW) computes one of 4 directional co-occurence matrices and features.
-    In debug=1 will return the 4 joint histograms for each ROI. This is used to compute texture
-    features for a pre-segmented region. The seg_co_occurence() method will be used for texture-
-    based segmentation.
-
-    Parameters 
-    ----------
-
-    label_image : {nd_array}
-        an image with labeled regions from get_blobs() method
-
-    raw_image : {nd_array}
-        raw image from which texture features get extracted 
-
-    distance : {int}
-        integer value of pixel offset in forming joint histogram. default value 2
-
-    orientation : {45, 90, 135, 180}
-        direction for pixel offet.
-
-    ROI : {dictionary}
-        Region of Interest structure that has blob bounding boxes. The largest
-        2D target bounding box is extracted.
-
-
-    Returns 
-    ----------
-
-    co_occurence_images : {dictionary}
-        contains 4 joint histogram images for each ROI 
-        returned if verbose=1
-
-    """
-
-    if orientation != 45 and orientation != 90 and orientation != 135 and orientation != 180: 
-        orientation = 90
-
-    epsilon  = 2.2e-16 
-    num_bits = 256
-    copy_image = raw_image.copy()
-    number_regions = ROI.size
-    indices = range(0, number_regions)
-    co_occurence_image_list = {}
-    for i in indices:
-        left   = ROI[i]['Left']
-        right  = ROI[i]['Right']
-        bottom = ROI[i]['Bottom']
-        top    = ROI[i]['Top']
-        Label  = ROI[i]['Label']
-        rows   = top-bottom
-        cols   = right-left
-        # copy the mask to section image
-        section = NP.zeros(rows*cols, dtype=label_image.dtype).reshape(rows, cols)
-        section[0:rows, 0:cols][label_image[bottom:top, left:right]==Label] = 1
-        source_region = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
-        cocm_block = NP.zeros(num_bits*num_bits, dtype=NP.int32).reshape(num_bits, num_bits)
-        source_region[0:rows, 0:cols] = copy_image[bottom:top, left:right] 
-        # scale segment to 8 bits. this needs to be smarter (e.g. use integrated histogram method)
-        max_value = source_region.max()
-        min_value = source_region.min()
-        scale = 255.0 / (max_value-min_value)
-        image_roi = (scale*(source_region-min_value)).astype(NP.int16)
-        # image_roi is short type
-        S.roi_co_occurence(section, image_roi, cocm_block, distance, orientation)
-        co_occurence_image_list[i] = cocm_block
-        # normalize the joint histogram prior to feature extraction
-        joint_histogram  = cocm_block.astype(NP.float64) 
-        joint_histogram  = joint_histogram / joint_histogram.sum()
-        # to prevent log(0)
-        joint_histogram += epsilon
-        # compute the com features
-        energy = joint_histogram.std()
-        H = joint_histogram * NP.log(joint_histogram)
-        entropy = H.sum()
-        r, c = joint_histogram.shape
-        [a, b] = NP.mgrid[1:c+1, 1:r+1]
-        contrast = ((NP.square(a-b))*joint_histogram).sum()
-        d = 1.0 + NP.abs(a-b)
-        homogeneity = (joint_histogram / d).sum()
-        ROI[i]['COM'][0] = distance
-        ROI[i]['COM'][1] = orientation 
-        ROI[i]['COM'][2] = energy
-        ROI[i]['COM'][3] = entropy
-        ROI[i]['COM'][4] = contrast
-        ROI[i]['COM'][5] = homogeneity
-
-    if verbose == 1:
-        return co_occurence_image_list
-    else:
-        return
-
-
-
-def region_grow(label_image, raw_image, ROI, roi_index, roi_inflate,
-                low_thresh=0.5, high_thresh=1.5, N_connectivity=3, debug=0):
-    """
-    region_grow(label_image, raw_image, ROI, roi_index, roi_inflate, stop_thresh)
-
-    starting from a seed (masks or regions in the label_image) grow the regions based
-    on (1) connectivity (2D or 3D) and (2) raw voxel values > stop threshold criterion.
-
-    Parameters 
-    ----------
-
-    label_image : {nd_array}
-        an image with labeled regions from get_blobs() method
-
-    raw_image : {nd_array}
-        raw image from which texture features get extracted 
-
-    ROI : {dictionary}
-        Region of Interest structure that has blob bounding boxes. The largest
-        2D target bounding box is extracted.
-
-    roi_index : {int}
-        the single ROI element to apply region growing to.
-
-    roi_inflate : {list}
-        the maximum increase in the ROI bounding box. For 3D the tuple is [layers, rows, cols]
-        and for 2D it is [rows, cols].
-
-    low_thresh : {float}
-        this is the percent of the voxel mean that the growing region must be GREATER than.
-        region growing terminates when the raw_image is BELOW this value.
-
-    high_thresh : {float}
-        this is the percent of the voxel mean that the growing region must be LESS than.
-        region growing terminates when the raw_image is ABOVE this value.
-
-    N_connectivity : {int}
-        for growing this indicates how connected in a 3x3 or 3x3x3 window the un-labeled
-        sample is. Make less than full connected for growing boundaries
-
-    Returns 
-    ----------
-
-    label : {nd_array}
-        the label image with the selected ROI after region growing. only returned
-        in debug mode.
-
-    """
-
-    _c_ext_struct = NP.dtype([('Left', 'i'),
-                              ('Right', 'i'),
-                              ('Top', 'i'),
-                              ('Bottom', 'i'),
-                              ('Front', 'i'),
-                              ('Back', 'i'),
-                              ('Label', 'i'),
-                              ('Mass', 'i'),
-                              ('cX', 'f'),
-                              ('cY', 'f'),
-                              ('cZ', 'f')]
-                             )
-
-    expanded_ROI = NP.zeros(1, dtype=_c_ext_struct)
-
-    dimensions  = label_image.ndim
-
-    if dimensions == 3:
-        z_ext = roi_inflate[0]
-        y_ext = roi_inflate[1]
-        x_ext = roi_inflate[2]
-        [layers, rows, cols]  = label_image.shape
-    else:
-        y_ext = roi_inflate[0]
-        x_ext = roi_inflate[1]
-        [rows, cols]  = label_image.shape
-
-    if dimensions == 2:  
-        left    = ROI[roi_index]['Left']-x_ext
-        right   = ROI[roi_index]['Right']+x_ext
-        bottom  = ROI[roi_index]['Bottom']-y_ext
-        top     = ROI[roi_index]['Top']+y_ext
-        Label   = ROI[roi_index]['Label']
-        lcutoff = low_thresh  * ROI[roi_index]['voxelMean']
-        hcutoff = high_thresh * ROI[roi_index]['voxelMean']
-        if left < 0: 
-            left = 0
-        if bottom < 0: 
-            bottom = 0
-        if right > cols-1: 
-            right = cols-1
-        if top > rows-1: 
-            top = rows-1
-        expanded_ROI['Left']   = left 
-        expanded_ROI['Right']  = right 
-        expanded_ROI['Top']    = top 
-        expanded_ROI['Bottom'] = bottom 
-        expanded_ROI['Label']  = Label 
-        rows    = top-bottom
-        cols    = right-left
-        label   = NP.zeros(rows*cols, dtype=NP.int16).reshape(rows, cols)
-        section = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
-        label   = label_image[bottom:top, left:right].copy()
-        section = (raw_image[bottom:top, left:right].astype(NP.float64)).copy()
-    elif dimensions == 3:  
-        left    = ROI[roi_index]['Left']-x_ext
-        right   = ROI[roi_index]['Right']+x_ext
-        bottom  = ROI[roi_index]['Bottom']-y_ext
-        top     = ROI[roi_index]['Top']+y_ext
-        front   = ROI[roi_index]['Front']-z_ext
-        back    = ROI[roi_index]['Back']+z_ext
-        Label   = ROI[roi_index]['Label']
-        lcutoff = low_thresh  * ROI[roi_index]['voxelMean']
-        hcutoff = high_thresh * ROI[roi_index]['voxelMean']
-        if left < 0: 
-            left = 0
-        if bottom < 0: 
-            bottom = 0
-        if right > cols-1: 
-            right = cols-1
-        if top > rows-1: 
-            top = rows-1
-        if front < 0: 
-            front = 0
-        if back > layers-1: 
-            back = layers-1
-        expanded_ROI['Left']   = left 
-        expanded_ROI['Right']  = right 
-        expanded_ROI['Top']    = top 
-        expanded_ROI['Bottom'] = bottom 
-        expanded_ROI['Back']   = back 
-        expanded_ROI['Front']  = front 
-        expanded_ROI['Label']  = Label 
-        rows    = top-bottom
-        cols    = right-left
-        layers  = back-front
-        label   = NP.zeros(layers*rows*cols, dtype=NP.int16).reshape(layers, rows, cols)
-        label   = label_image[front:back, bottom:top, left:right].copy()
-        section = NP.zeros(layers*rows*cols, dtype=NP.float64).reshape(layers, rows, cols)
-        section = (raw_image[front:back, bottom:top, left:right].astype(NP.float64)).copy()
-
-    #
-    # this newgrow_ROI gets filled in and the label image is grown
-    #
-
-    newgrow_ROI = NP.zeros(1, dtype=_c_ext_struct)
-    S.region_grow(section, label, expanded_ROI, newgrow_ROI, lcutoff, hcutoff, Label, N_connectivity)
-
-    if debug==1:  
-        #
-        # do not update ROI for index and the label_image 
-        #
-        return label
-
-    else:
-        #
-        # update (overwrite) ROI for index and the label_image 
-        #
-        if dimensions == 2:  
-            ROI[roi_index]['Left']   = newgrow_ROI['Left']
-            ROI[roi_index]['Right']  = newgrow_ROI['Right']
-            ROI[roi_index]['Top']    = newgrow_ROI['Top']
-            ROI[roi_index]['Bottom'] = newgrow_ROI['Bottom']
-            left   = ROI[roi_index]['Left']
-            right  = ROI[roi_index]['Right']
-            top    = ROI[roi_index]['Top']
-            bottom = ROI[roi_index]['Bottom']
-            rows   = top-bottom
-            cols   = right-left
-            label_image[bottom:top,left:right] = label[0:rows,0:cols]
-        elif dimensions == 3:  
-            ROI[roi_index]['Left']   = newgrow_ROI['Left']
-            ROI[roi_index]['Right']  = newgrow_ROI['Right']
-            ROI[roi_index]['Top']    = newgrow_ROI['Top']
-            ROI[roi_index]['Bottom'] = newgrow_ROI['Bottom']
-            ROI[roi_index]['Front']  = newgrow_ROI['Front']
-            ROI[roi_index]['Back']   = newgrow_ROI['Back']
-            left   = expanded_ROI['Left']
-            right  = expanded_ROI['Right']
-            top    = expanded_ROI['Top']
-            bottom = expanded_ROI['Bottom']
-            front  = expanded_ROI['Front']
-            back   = expanded_ROI['Back']
-            rows   = top-bottom
-            cols   = right-left
-            layers = back-front
-            label_image[front:back,bottom:top,left:right] = label[0:layers,0:rows,0:cols]
-                             
-        return 
-
-
-
-def seg_co_occurence(raw_image, window=16, distance=2, orientation=90):
-    """
-    cocm_images = seg_co_occurence(raw_image, window=16, distance=2, orientation=90)
-
-    (N-S, E-W, NW-SE, NE-SW) computes one of 4 directional co-occurence matrices and features.
-    In debug=1 will return the 4 joint histograms for each ROI.
-
-    The seg_co_occurence() method is used for texture-based segmentation. Feature images are
-    returned from which segmentation can be later performed.
-
-    ****
-    NOTE: This is very slow and a fast method using Unsers histogram approximation will be 
-    added in the future.
-    ****
-
-
-    Parameters 
-    ----------
-
-    raw_image : {nd_array}
-        raw image from which texture features get extracted 
-
-    window : {int}
-        integer value of moving 2D window. Window slides in 2D over image and is the
-        region-of-interest from which co-occurence texture features are extracted. The
-        window is 2D square so only a single value is entered. Default window is 32x32. 
-
-    distance : {int}
-        integer value of pixel offset in forming joint histogram. default value 2
-
-    orientation : {45, 90, 135, 180}
-        direction for pixel offet.
-
-    Returns 
-    ----------
-
-    cocm_images : {dictionary}
-        
-        co_occurence_feature_images. contains 4 normalized feature
-        windows with keys: energy, entropy, contrast and homogeneity.
-
-    """
-
-    if orientation != 45 and orientation != 90 and orientation != 135 and orientation != 180: 
-        orientation = 90
-
-    epsilon      = 2.2e-16 
-    num_bits     = 256
-    copy_image   = raw_image.copy()
-    [rows, cols] = copy_image.shape
-    row_indices  = range(window, rows-window)
-    col_indices  = range(window, cols-window)
-
-    # create a fixed mask and scratch window for raw source
-    section       = NP.ones(2*window*2*window, dtype=NP.int16).reshape(2*window, 2*window)
-    source_region = NP.zeros(2*window*2*window, dtype=NP.float64).reshape(2*window, 2*window)
-
-    # output images
-    energy_image      = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
-    entropy_image     = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
-    homogeneity_image = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
-    contrast_image    = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
-    cocm_block        = NP.zeros(num_bits*num_bits, dtype=NP.int32).reshape(num_bits, num_bits)
-    
-    for i in row_indices:
-        bottom = i - window
-        top    = i + window
-        for j in col_indices:
-            left  = j - window
-            right = j + window 
-            source_region[0:2*window, 0:2*window] = copy_image[bottom:top, left:right] 
-            # scale segment to 8 bits. this needs to be smarter (e.g. use integrated histogram method)
-            max_value = source_region.max()
-            min_value = source_region.min()
-            scale     = 255.0 / (max_value-min_value)
-            image_roi = (scale*(source_region-min_value)).astype(NP.int16)
-            # image_roi is short type
-            cocm_block[:] = 0.0
-            S.roi_co_occurence(section, image_roi, cocm_block, distance, orientation)
-            # normalize the joint histogram prior to feature extraction
-            joint_histogram = cocm_block.astype(NP.float64) 
-            joint_histogram = joint_histogram / joint_histogram.sum()
-            # to prevent log(0)
-            joint_histogram += epsilon
-            # compute the com features
-            energy      = joint_histogram.std()
-            H           = joint_histogram * NP.log(joint_histogram)
-            entropy     = H.sum()
-            r, c        = joint_histogram.shape
-            [a, b]      = NP.mgrid[1:c+1, 1:r+1]
-            contrast    = ((NP.square(a-b))*joint_histogram).sum()
-            d           = 1.0 + NP.abs(a-b)
-            homogeneity = (joint_histogram / d).sum()
-            # store the feature pixel for the 4 images
-            energy_image[i, j]      = energy
-            entropy_image[i, j]     = entropy
-            contrast_image[i, j]    = contrast
-            homogeneity_image[i, j] = homogeneity
-
-    scale_energy      = 1.0 / max(energy_image.max(), abs(energy_image.min()))
-    scale_entropy     = 1.0 / max(entropy_image.max(), abs(entropy_image.min()))
-    scale_contrast    = 1.0 / max(contrast_image.max(), abs(contrast_image.min()))
-    scale_homogeneity = 1.0 / max(homogeneity_image.max(), abs(homogeneity_image.min()))
-
-    energy_image      = scale_energy      * energy_image
-    entropy_image     = scale_entropy     * entropy_image
-    homogeneity_image = scale_homogeneity * homogeneity_image
-    contrast_image    = scale_contrast    * contrast_image
-
-    cocm_images = {'energy_image' : energy_image,  'entropy_image' : entropy_image, 
-                   'homogeneity_image' : homogeneity_image,  'contrast_image' : contrast_image} 
-
-    return cocm_images 
-
-
-def roi_mat_filter(label_image, thin_kernel, ROI):
-    """
-    thin_edge_image = roi_mat_filter(label_image, thin_kernel, ROI)
-
-    gets the largest object in the ROI list and returns the medial axis for
-    that object. Idea is that the largest object is a reference, e.g. the skull
-    for anatomical MRI.
-
-    Parameters 
-    ----------
-
-    label_image : {nd_array}
-        an image with labeled regions from get_blobs() method
-
-    thin_kernel : {dictionary}
-        set of 8 'J' and 'K' 3x3 masks from build_morpho_thin_masks() method
-
-    ROI : {dictionary}
-        Region of Interest structure that has blob bounding boxes. The largest
-        2D target bounding box is extracted.
-
-    Returns 
-    ----------
-
-    thin_edge_image : {nd_array}
-        thinned edge image for the largest object.
-
-    """
-
-
-    [rows, cols] = label_image.shape
-    # destination image
-    thin_edge_image = NP.zeros(rows*cols, dtype=NP.uint16).reshape(rows, cols)
-    # scratch memory for thin 
-    input           = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
-    cinput          = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
-    erosion         = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
-    dialation       = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
-    hmt             = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
-    copy            = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
-
-    bbox = get_max_bounding_box(ROI)
-
-    left     = bbox['Left']-1
-    right    = bbox['Right']+1
-    bottom   = bbox['Bottom']-1
-    top      = bbox['Top']+1
-    Label    = bbox['Label']
-
-    if left < 0: 
-        left = 0
-    if bottom < 0: 
-        bottom = 0
-    if right > cols-1: 
-        right = cols-1
-    if top > rows-1: 
-        top = rows-1
-
-    inflate  = 1
-    roi_rows = top-bottom+2*inflate
-    roi_cols = right-left+2*inflate
-    rgrows   = top-bottom
-    rgcols   = right-left
-    # clear the memory
-    input[0:roi_rows, 0:roi_cols] = 0
-    # load the labeled region 
-    input[inflate:inflate+rgrows, inflate:inflate+rgcols] \
-         [label_image[bottom:top, left:right]==Label] = 1 
-    # thin this region
-    S.thin_filter(thin_kernel['jmask'], thin_kernel['kmask'], thin_kernel['number3x3Masks'],
-                  roi_rows, roi_cols, cols, input, cinput, erosion, dialation, hmt, copy)
-
-    # accumulate the images (do not over-write). for overlapping regions
-    input[inflate:rgrows+inflate,inflate:rgcols+inflate] \
-         [input[inflate:rgrows+inflate,inflate:rgcols+inflate]==1] = Label 
-    thin_edge_image[bottom:top,left:right] = input[inflate:rgrows+inflate,inflate:rgcols+inflate] 
-
-    return thin_edge_image
-
-def mat_filter(label_image, thin_kernel, ROI=None):
-    """
-    mat_image = mat_filter(label_image, thin_kernel, ROI=None)
-
-    takes the ROI dictionary with the blob bounding boxes and thins each blob
-    giving the medial axis. if ROI is null will create a single ROI with the
-    bounding box set equal to the full image
-
-    Parameters 
-    ----------
-
-    label_image : {nd_array}
-        an image with labeled regions from get_blobs() method
-
-    thin_kernel : {dictionary}
-        set of 8 'J' and 'K' 3x3 masks from build_morpho_thin_masks() method
-
-    ROI : {dictionary}
-        Region of Interest structure that has blob bounding boxes
-
-    Returns 
-    ----------
-
-    mat_image : {nd_array}
-        thinned edge image
-
-    """
-    if ROI==None:
-        ROIList = NP.zeros(1, dtype=_objstruct)
-        [rows, cols] = label_image.shape
-        ROIList['Left']   = 2
-        ROIList['Right']  = cols-3
-        ROIList['Bottom'] = 2
-        ROIList['Top']    = rows-3
-
-    [rows, cols] = label_image.shape
-    # destination image
-    thin_edge_image = NP.zeros(rows*cols, dtype=NP.uint16).reshape(rows, cols)
-    mat_image       = NP.zeros(rows*cols, dtype=NP.uint16).reshape(rows, cols)
-    # scratch memory for thin 
-    input           = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
-    cinput          = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
-    erosion         = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
-    dialation       = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
-    hmt             = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
-    copy            = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
-
-    number_regions = ROI.size
-    indices = range(0, number_regions)
-    inflate = 1
-    for i in indices:
-        left     = ROI[i]['Left']-1
-        right    = ROI[i]['Right']+1
-        bottom   = ROI[i]['Bottom']-1
-        top      = ROI[i]['Top']+1
-        Label    = ROI[i]['Label']
-        if left < 0: 
-            left = 0
-        if bottom < 0: 
-            bottom = 0
-        if right > cols-1: 
-            right = cols-1
-        if top > rows-1: 
-            top = rows-1
-
-        roi_rows = top-bottom+2*inflate
-        roi_cols = right-left+2*inflate
-        rgrows   = top-bottom
-        rgcols   = right-left
-        # clear the memory
-        input[0:roi_rows, 0:roi_cols] = 0
-        # load the labeled region 
-        input[inflate:inflate+rgrows, inflate:inflate+rgcols] \
-             [label_image[bottom:top, left:right]==Label] = 1 
-        # thin this region
-        S.thin_filter(thin_kernel['jmask'], thin_kernel['kmask'], thin_kernel['number3x3Masks'],
-                      roi_rows, roi_cols, cols, input, cinput, erosion, dialation, hmt, copy)
-
-        # accumulate the images (do not over-write). for overlapping regions
-        input[inflate:rgrows+inflate,inflate:rgcols+inflate] \
-             [input[inflate:rgrows+inflate,inflate:rgcols+inflate]==1] = Label 
-        thin_edge_image[bottom:top,left:right] = thin_edge_image[bottom:top,left:right] + \
-                                                 input[inflate:rgrows+inflate,inflate:rgcols+inflate] 
-            
-
-    # accumulate overlaps set back to binary at later date
-    mat_image[:, :] = thin_edge_image[:, :]
-
-    return mat_image
-
-
-def laws_texture_filter(raw_image, label_image, laws_kernel, ROI=None, dc_thres=1.0,
-                        mean_feature=1, verbose=0):
-    """
-    texture_images = laws_texture_filter(raw_image, label_image, laws_kernel, ROI=None, verbose=1)
-    .
-        OR
-    .
-    laws_texture_filter(raw_image, label_image, laws_kernel, ROI=None, verbose=0)
-
-    Parameters 
-    ----------
-
-    raw_image : {nd_array}
-        raw double image 
-
-    label_image : {nd_array}
-        an image with labeled regions from get_blobs() method
-
-    laws_kernel : {dictionary}
-        set of 6 length-7 Law's texture feature kernels 
-
-    ROI : {dictionary}
-        Region of Interest structure that has blob bounding boxes
-
-    dc_thres : {float}
-        used as a filter. Sets texture feature to 0.0 when the 
-        mean level is above this. Removes the low frequency, high amplitude
-        image regions from the feature list
-
-    mean_feature : {0, 1}, optional
-        when set to 1, the feature is the mean value of the
-        selected Law's texture filter. When 0 the feature is
-        the standard deviation.
-
-    verbose : {0, 1}, optional
-        determines if return is to include Law's filter images
-
-    Returns 
-    ----------
-
-    laws_image : {dictionary}
-        contains 21 Laws filtered  regions for each ROI 
-        returned if verbose=1
-        
-
-    """
-    if ROI==None:
-        ROI= NP.zeros(1, dtype=_objstruct)
-        [rows, cols] = label_image.shape
-        ROI['Left']   = 2
-        ROI['Right']  = cols-3
-        ROI['Bottom'] = 2
-        ROI['Top']    = rows-3
-
-    laws_image_list = {}
-    number_regions  = ROI.size
-    layers          = laws_kernel['filters']
-    indices         = range(0, number_regions)
-    filters         = range(0, layers)
-    for i in indices:
-        left   = ROI[i]['Left']
-        right  = ROI[i]['Right']
-        bottom = ROI[i]['Bottom']
-        top    = ROI[i]['Top']
-        Label  = ROI[i]['Label']
-        rows   = top-bottom
-        cols   = right-left
-        label_region  = NP.zeros(rows*cols, dtype=NP.uint16).reshape(rows, cols)
-        source_region = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
-        laws_block    = NP.zeros(layers*rows*cols, dtype=NP.float32).reshape(layers, rows, cols)
-        # load the labeled region 
-        label_region[0:rows,  0:cols][label_image[bottom:top, left:right]==Label] = 1 
-        source_region[0:rows, 0:cols] = raw_image[bottom:top, left:right] 
-
-        S.laws_texture_metric(label_region, source_region, laws_block, laws_kernel['numKernels'],
-                              laws_kernel['kernelSize'], laws_kernel['filters'],
-                              laws_kernel['coefficients'][0], laws_kernel['coefficients'][1],
-                              laws_kernel['coefficients'][2], laws_kernel['coefficients'][3],
-                              laws_kernel['coefficients'][4], laws_kernel['coefficients'][5])
-
-        for j in filters:
-            # compute the energy measure for each filter in the ROI
-            mask_image = laws_block[j, :, :][label_region[:, :]>0]
-            mean = abs(mask_image.mean())
-            std  = mask_image.std()
-            if mean > dc_thres:
-                mean = 0.0
-                std = 0.0
-            if mean_feature == 1:
-                ROI[i]['TEM'][j] = mean 
-            else:
-                ROI[i]['TEM'][j] = std 
-
-        ROI[i]['TEM'][:] = ROI[i]['TEM'][:] / ROI[i]['TEM'][:].max() 
-        # accumulate the 21 Law's filtered ROI's and optional
-        # return as image (3D)
-        laws_image_list[i] = laws_block
-
-    if verbose == 1:
-        return laws_image_list
-    else:
-        return 
-
-
-
-def get_voxel_measures(label_image, raw_image, ROI=None):
-    """
-    mat_image = mat_filter(label_image, raw_image, ROI=None)
-
-    takes the ROI dictionary with the blob bounding boxes and gets the voxel measures
-    from each ROI in the raw data.
-
-    Parameters 
-    ----------
-
-    label_image : {nd_array}
-        an image with labeled regions from get_blobs() method
-
-    raw_image : {nd_array}
-        the original double image (raw voxels) from which voxel measures are made
-
-    ROI : {dictionary}
-        Region of Interest structure that has blob bounding boxes
-
-
-    Returns 
-    ----------
-
-    none
-
-    """
-
-    dimensions = label_image.ndim
-
-    if ROI==None:
-        ROIList = NP.zeros(1, dtype=_objstruct)
-        if dimensions == 2:  
-            [rows, cols] = label_image.shape
-            ROIList['Left']   = 1
-            ROIList['Right']  = cols-1
-            ROIList['Bottom'] = 1
-            ROIList['Top']    = rows-1
-        elif dimensions == 3:  
-            [layers, rows, cols] = label_image.shape
-            ROIList['Left']   = 1
-            ROIList['Right']  = cols-1
-            ROIList['Bottom'] = 1
-            ROIList['Top']    = rows-1
-            ROIList['Front']  = 1
-            ROIList['Back']   = layers-1
-
-    number_regions = ROI.size
-    indices = range(0, number_regions)
-    inflate = 1
-    for i in indices:
-        if dimensions == 2:  
-            left   = ROI[i]['Left']
-            right  = ROI[i]['Right']
-            bottom = ROI[i]['Bottom']
-            top    = ROI[i]['Top']
-            Label  = ROI[i]['Label']
-            rows   = top-bottom-1
-            cols   = right-left-1
-            section= NP.zeros(rows*cols, dtype=raw_image.dtype).reshape(rows, cols)
-            section = raw_image[bottom:top, left:right] \
-                               [label_image[bottom:top, left:right]==Label]
-        elif dimensions == 3:  
-            left   = ROI[i]['Left']
-            right  = ROI[i]['Right']
-            bottom = ROI[i]['Bottom']
-            top    = ROI[i]['Top']
-            front  = ROI[i]['Front']
-            back   = ROI[i]['Back']
-            Label  = ROI[i]['Label']
-            rows   = top-bottom-1
-            cols   = right-left-1
-            layers = back-front-1
-            section= NP.zeros(layers*rows*cols, dtype=raw_image.dtype).reshape(layers, rows, cols)
-            section = raw_image[front:back, bottom:top, left:right] \
-                               [label_image[front:back, bottom:top, left:right]==Label]
-
-        mask = section[section>0]
-        ROI[i]['voxelMean'] = mask.mean()
-        ROI[i]['voxelVar']  = mask.std()
-
-    return 
-
-
-def get_blob_regions(labeled_image, groups, dust=16):
-    """
-    ROIList = get_blob_regions(labeled_image, groups, dust=16)
-
-    get the bounding box for each labelled blob in the image
-    allocate the dictionary structure that gets filled in with later
-    stage processing to add blob features.
-
-    Parameters 
-    ----------
-
-    label_image : {nd_array}
-        a 2D or 3D image with labeled regions from get_blobs() method
-
-    groups : {int}
-        number of blobs in image determined by get_blobs() method
-
-    Returns 
-    ----------
-
-    ROIList : {dictionary}
-        structure that has the bounding box and area of each blob
-
-
-
-    """
-
-    _c_ext_struct = NP.dtype([('Left', 'i'),
-                              ('Right', 'i'),
-                              ('Top', 'i'),
-                              ('Bottom', 'i'),
-                              ('Front', 'i'),
-                              ('Back', 'i'),
-                              ('Label', 'i'),
-                              ('Mass', 'i'),
-                              ('cX', 'f'),
-                              ('cY', 'f'),
-                              ('cZ', 'f')]
-                             )
-
-    c_ext_ROI = NP.zeros(groups, dtype=_c_ext_struct)
-    ROIList = NP.zeros(groups, dtype=_objstruct)
-    # return the bounding box for each connected edge
-    S.get_blob_regions(labeled_image, c_ext_ROI)
-
-    indices = range(0, groups)
-    for i in indices:
-        ROIList[i]['Left']   = c_ext_ROI[i]['Left']
-        ROIList[i]['Right']  = c_ext_ROI[i]['Right']
-        ROIList[i]['Bottom'] = c_ext_ROI[i]['Bottom']
-        ROIList[i]['Top']    = c_ext_ROI[i]['Top']
-        ROIList[i]['Front']  = c_ext_ROI[i]['Front']
-        ROIList[i]['Back']   = c_ext_ROI[i]['Back']
-        ROIList[i]['Label']  = c_ext_ROI[i]['Label']
-        ROIList[i]['Mass']   = c_ext_ROI[i]['Mass']
-        ROIList[i]['cX']     = c_ext_ROI[i]['cX']
-        ROIList[i]['cY']     = c_ext_ROI[i]['cY']
-        ROIList[i]['cZ']     = c_ext_ROI[i]['cZ']
-
-    return ROIList[ROIList['Mass']>dust]
-
-
-def get_blobs(binary_edge_image, mask=1):
-    """
-
-    labeled_edge_image, groups = get_blobs(binary_edge_image)
-
-    get the total number of blobs in a 2D or 3D image and convert the 
-    binary image (or volume) to labelled regions
-
-    Parameters 
-    ----------
-
-    binary_edge_image : {nd_array}
-        an binary image/volume
-
-    mask : {int}
-        the size of the 2D or 3D connectivity mask. For 2D this is 1, 4 or 8.
-        For 3D this is 1, 6, 14 or 28. Mask = 1 is ANY connection in 3x3
-        or 3x3x3 mask for 2D or 3D, respectively.
-
-    Returns 
-    ----------
-
-    label_image : {nd_array}
-        an image/volume with labeled regions from get_blobs() method
-
-    groups : {int}
-        number of blobs in image determined by get_blobs() method
-
-    """
-
-    dimensions = binary_edge_image.ndim
-    if dimensions == 2:  
-        if mask != 1 and mask != 4 and mask != 8:
-            mask = 1 
-        [rows, cols] = binary_edge_image.shape
-        labeled_edge_image_or_vol = NP.zeros(rows*cols, dtype=NP.uint16).reshape(rows, cols)
-    elif dimensions == 3:
-        if mask != 1 and mask != 6 and mask != 14 and mask != 28:
-            mask = 1 
-        [layers, rows, cols] = binary_edge_image.shape
-        labeled_edge_image_or_vol = NP.zeros(layers*rows*cols, dtype=NP.uint16).reshape(layers, rows, cols)
-    else:
-        labeled_edge_image_or_vol = None
-        groups = 0
-        return labeled_edge_image_or_vol, groups
-
-    groups = S.get_blobs(binary_edge_image, labeled_edge_image_or_vol, mask)
-
-    return labeled_edge_image_or_vol, groups
-
-
-def sobel_edges(sobel_edge_image, sobel_stats, mode=1, sobel_threshold=0.3):
-    """
-    sobel_edge = sobel_edges(sobel_edge_image, sobel_stats, mode=1, sobel_threshold=0.3)
-    take sobel-filtered image and return binary edges
-
-    Parameters 
-    ----------
-
-    sobel_edge_image : {nd_array}
-        edge-filtered image from sobel_image() method
-
-    sobel_stats : {dictionary}
-        mean and nonzero min, max of sobel filtering
-
-    mode : {0, 1}, optional
-        threshold based on histogram mean(0) or mode(1)
-
-    sobel_threshold : {float}, optional
-        low threshold applied to edge filtered image for edge generation
-
-    Returns 
-    ----------
-
-    sobel_edge : {nd_array}
-        binary edge-image
-
-    """
-    [rows, cols] = sobel_edge_image.shape
-    sobel_edge = NP.zeros(rows*cols, dtype=NP.uint16).reshape(rows, cols)
-    S.sobel_edges(sobel_edge_image, sobel_edge, sobel_stats['ave_gt0'], sobel_stats['min_gt0'],
-                  sobel_stats['max_gt0'], mode, sobel_threshold)
-
-    return sobel_edge
-
-
-def sobel_image(filtered_slice):
-    """
-    sobel_edge_image, sobel_stats = sobel_image(filtered_slice)
-
-    take 2D raw or filtered image slice and get sobel-filtered image 
-
-    Parameters 
-    ----------
-
-    filtered_slice : {nd_array}
-        raw or pre-processed (filtered and thresholded) 2D image 
-
-    Returns 
-    ----------
-
-    sobel_edge_image : {nd_array}
-        edge-filtered image from sobel_image() method
-
-    sobel_stats : {dictionary}
-        mean and nonzero min, max of sobel filtering
-
-    """
-    filtered_slice = filtered_slice.astype(NP.float64)
-    [rows, cols] = filtered_slice.shape
-    sobel_edge_image = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
-    pAve, min_value, max_value = S.sobel_image(filtered_slice, sobel_edge_image)
-
-    # can replace this with numpy calls for the image stats. but C ext is faster
-    #   S.sobel_image(filtered_slice, sobel_edge_image)
-    #   sobel_mask = sobel_edge_image[sobel_edge_image>0]
-    #   pAve = sobel_mask.mean()
-    #   min_value = sobel_mask.min()
-    #   max_value = sobel_mask.max()
-
-    sobel_stats= {'ave_gt0' : pAve, 'min_gt0': min_value, 'max_gt0': max_value} 
-
-    return sobel_edge_image, sobel_stats
-
-def pre_filter(slice, filter, low_threshold=0, high_threshold=0, conv_binary=0):
-    """
-    edge_filter = pre_filter(slice, filter, low_threshold=0, high_threshold=slice.max(), conv_binary=0)
-
-    take 2D image slice and filter and pre-filter and threshold prior to segmentation
-
-    Parameters 
-    ----------
-    slice : {nd_array}
-        input 2D image. gets cast to int16
-
-    filter : {dictionary}
-        2D filter kernel set from build_2d_kernel()
-
-    low_threshold : {int}, optional
-        default 0
-
-    high_threshold : {int}, optional
-        default max value of source image
-
-    conv_binary : {0, 1}, optional
-        flag to convert edge_filter image to binary valued. default 
-        is binary conversion off
-
-    Returns 
-    ----------
-    edge_filter : {nd_array}
-        filtered and thresholded image that can be (optional) binary.
-
-    """
-
-    # make sure the input is 16 bits. this is input to edge machine
-    # so can handle raw and 8 bit scaled inputs
-    if high_threshold==0:
-        # default to the maximum value of the image
-        high_threshold = slice.max()
-
-    slice = slice.astype(NP.int16)
-    [rows, cols] = slice.shape
-    edge_image = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
-    S.edge_prefilter(low_threshold, high_threshold, filter['kernelSize'], filter['kernel'],
-                     slice, edge_image)
-
-    if conv_binary == 1:
-        edge_image[edge_image>0] = 1
-        edge_image = edge_image.astype(NP.uint16)
-
-    return edge_image
-
-def get_max_bounding_box(ROI):
-    """
-    bounding_box = get_max_bounding_box(ROI)
-
-    take an ROI structure and find the maximum area bounding box
-
-    Parameters 
-    ----------
-
-    ROI : {dictionary}
-        the ROI is the automatically extracted blob regions of interest
-        and contains the rectangular bounding box of each blob.
-
-    Returns 
-    ----------
-
-    bounding_box : {dictionary}
-        the Left, Right, Top Bottom and Label of the LARGEST bounding box in the ROI
-
-    """
-    max_index = ROI[:]['Mass'].argmax()
-    bounding_box = {'Left' : ROI[max_index]['Left'], 'Right' : ROI[max_index]['Right'],
-                    'Top' : ROI[max_index]['Top'], 'Bottom' : ROI[max_index]['Bottom'],
-                    'Label' : ROI[max_index]['Label']} 
-
-    return bounding_box 
-
-def get_all_bounding_boxes(ROI):
-    """
-    measures = get_all_bounding_boxes(ROI)
-
-    get all bounding boxes in the ROI (feature) dictionary
-
-    Parameters 
-    ----------
-
-    ROI : {dictionary}
-        the ROI is the automatically extracted blob regions of interest
-        and contains the rectangular bounding box of each blob.
-
-    Returns 
-    ----------
-
-    measures : {dictionary}
-        the Left, Right, Top and Bottom of all bounding boxes in the ROI
-
-    """
-
-    number = ROI.size
-    indices = range(0, ROI.size)
-    _shortstruct = NP.dtype([('left', 'i'),
-                             ('right', 'i'),
-                             ('top', 'i'),
-                             ('bottom', 'i')])
-    measures = NP.zeros(number, dtype=_shortstruct)
-    for i in indices:
-        measures[i]['left']   = ROI[i]['Left']
-        measures[i]['right']  = ROI[i]['Right']
-        measures[i]['top']    = ROI[i]['Top']
-        measures[i]['bottom'] = ROI[i]['Bottom']
-
-    return measures
-
-
-def build_2d_kernel(aperature=21, hiFilterCutoff=10.0):
-    """
-
-    FIRFilter = build_2d_kernel(aperature, hiFilterCutoff)
-
-    build hamming-windowed FIR filter with sinc kernel
-
-    Parameters 
-    ----------
-
-    aperature : {int}, optional
-        the number of coefficients in the filter. default is 21. needs to be ODD
-
-    hiFilterCutoff : {float}
-        the upper cutoff in digital frequency units
-
-    Returns 
-    ----------
-
-    FIRFilter : {dictionary}
-        filter kernel
-
-    """
-
-    rad = math.pi / 180.0
-    HalfFilterTaps = (aperature-1) / 2
-    kernel = NP.zeros((aperature), dtype=NP.float64)
-    LC = 0.0
-    HC = hiFilterCutoff * rad 
-    t2 = 2.0 * math.pi
-    t1 = 2.0 * HalfFilterTaps + 1.0
-    indices = range(-HalfFilterTaps, HalfFilterTaps+1, 1)  
-    j = 0
-    for i in indices:
-        if i == 0:
-            tLOW  = LC
-            tHIGH = HC
-        else:
-            tLOW  = math.sin(i*LC)/i
-            tHIGH = math.sin(i*HC)/i
-        # Hamming window
-        t3 = 0.54 + 0.46*(math.cos(i*t2/t1))
-        t4 = t3*(tHIGH-tLOW)
-        kernel[j] = t4
-        j += 1
-
-    # normalize the kernel
-    sum = kernel.sum()
-    kernel /= sum
-
-    FIRFilter= {'kernelSize' : HalfFilterTaps, 'kernel': kernel} 
-
-    return FIRFilter
-
-def build_d_gauss_kernel(gWidth=20, sigma=1.0):
-
-    """
-    DGFilter = build_d_gauss_kernel(gWidth, sigma)
-
-    build the derivative of Gaussian kernel for Canny edge filter
-
-    Parameters 
-    ----------
-    gWdith : {int}, optional
-         width of derivative of Gaussian kernel.
-         default value is 20
-
-    sigma : {float}, optional
-        sigma term of derivative of Gaussian kernel
-         default value is 1.0
-
-    Returns 
-    ----------
-
-    DGFilter : {dictionary}
-        filter kernel
-
-    """
-
-    kernel  = NP.zeros((1+gWidth), dtype=NP.float64)
-    indices = range(0, gWidth)  
-
-    for i in indices:
-        kernel[i]  = math.exp(((-i*i)/(2.0 * sigma * sigma)))
-        kernel[i] *= -(i / (sigma * sigma))
-
-    DGFilter= {'kernelSize' : gWidth, 'coefficients': kernel} 
-
-    return DGFilter
-
-def build_morpho_thin_masks():
-
-    """
-    MATFilter = build_morpho_thin_masks()
-
-    build 2 sets (J and K) of 8 3x3 morphology masks (structuring elements)
-    to implement thinning (medial axis transformation - MAT)
-
-
-    Parameters 
-    ----------
-
-    None
-
-    Returns 
-    ----------
-
-    MATFilter : {dictionary}
-        morphology filter kernels. there are 2 sets of 8 3x3 masks
-
-    """
-
-    # (layers, rows, cols)
-    size = (8*3*3)
-    J_mask = NP.zeros(size, dtype=NP.int16)
-    K_mask = NP.zeros(size, dtype=NP.int16)
-
-    maskCols = 3
-    # load the 8 J masks for medial axis transformation
-    Column = 0
-    J_mask[0+maskCols*(Column+0)] = 1
-    J_mask[0+maskCols*(Column+1)] = 1
-    J_mask[0+maskCols*(Column+2)] = 1
-    J_mask[1+maskCols*(Column+1)] = 1
-
-    Column += 3
-    J_mask[0+maskCols*(Column+1)] = 1
-    J_mask[1+maskCols*(Column+1)] = 1
-    J_mask[1+maskCols*(Column+2)] = 1
-
-    Column += 3
-    J_mask[0+maskCols*(Column+0)] = 1
-    J_mask[1+maskCols*(Column+0)] = 1
-    J_mask[2+maskCols*(Column+0)] = 1
-    J_mask[1+maskCols*(Column+1)] = 1
-
-    Column += 3
-    J_mask[0+maskCols*(Column+1)] = 1
-    J_mask[1+maskCols*(Column+0)] = 1
-    J_mask[1+maskCols*(Column+1)] = 1
-
-    Column += 3
-    J_mask[0+maskCols*(Column+2)] = 1
-    J_mask[1+maskCols*(Column+1)] = 1
-    J_mask[1+maskCols*(Column+2)] = 1
-    J_mask[2+maskCols*(Column+2)] = 1
-
-    Column += 3
-    J_mask[1+maskCols*(Column+0)] = 1
-    J_mask[1+maskCols*(Column+1)] = 1
-    J_mask[2+maskCols*(Column+1)] = 1
-
-    Column += 3
-    J_mask[1+maskCols*(Column+1)] = 1
-    J_mask[2+maskCols*(Column+0)] = 1
-    J_mask[2+maskCols*(Column+1)] = 1
-    J_mask[2+maskCols*(Column+2)] = 1
-
-    Column += 3
-    J_mask[1+maskCols*(Column+1)] = 1
-    J_mask[1+maskCols*(Column+2)] = 1
-    J_mask[2+maskCols*(Column+1)] = 1
-
-    # load the 8 K masks for medial axis transformation
-    Column = 0
-    K_mask[2+maskCols*(Column+0)] = 1
-    K_mask[2+maskCols*(Column+1)] = 1
-    K_mask[2+maskCols*(Column+2)] = 1
-    
-    Column += 3
-    K_mask[1+maskCols*(Column+0)] = 1
-    K_mask[2+maskCols*(Column+0)] = 1
-    K_mask[2+maskCols*(Column+1)] = 1
-    
-    Column += 3
-    K_mask[0+maskCols*(Column+2)] = 1
-    K_mask[1+maskCols*(Column+2)] = 1
-    K_mask[2+maskCols*(Column+2)] = 1
-
-    Column += 3
-    K_mask[1+maskCols*(Column+2)] = 1
-    K_mask[2+maskCols*(Column+1)] = 1
-    K_mask[2+maskCols*(Column+2)] = 1
-
-    Column += 3
-    K_mask[0+maskCols*(Column+0)] = 1
-    K_mask[1+maskCols*(Column+0)] = 1
-    K_mask[2+maskCols*(Column+0)] = 1
-
-    Column += 3
-    K_mask[0+maskCols*(Column+1)] = 1
-    K_mask[0+maskCols*(Column+2)] = 1
-    K_mask[1+maskCols*(Column+2)] = 1
-
-    Column += 3
-    K_mask[0+maskCols*(Column+0)] = 1
-    K_mask[0+maskCols*(Column+1)] = 1
-    K_mask[0+maskCols*(Column+2)] = 1
-
-    Column += 3
-    K_mask[0+maskCols*(Column+0)] = 1
-    K_mask[0+maskCols*(Column+1)] = 1
-    K_mask[1+maskCols*(Column+0)] = 1
-
-    MATFilter = {'number3x3Masks' : 8, 'jmask' : J_mask, 'kmask' : K_mask} 
-
-    return MATFilter 
-
-
-def build_laws_kernel():
-
-    """
-    LAWSFilter = build_laws_kernel()
-
-    build 6 length-7 Law's texture filter masks
-    mask names are: 'L', 'S', 'E', 'W', 'R', 'O'
-
-
-    Parameters 
-    ----------
-
-    None
-
-    Returns 
-    ----------
-
-    LAWSFilter : {dictionary}
-        a set of 6 length-7 Laws texture kernels
-
-    """
-    aperature = (6, 7)
-    coefficients = NP.zeros((aperature), dtype=NP.float64)
-    names = ('L', 'E', 'S', 'W', 'R', 'O' )
-
-    coefficients[0, :] =  ( 1.0,  6.0,  15.0, 20.0,  15.0,  6.0,  1.0)
-    coefficients[1, :] =  (-1.0, -4.0,  -5.0,  0.0,   5.0,  4.0,  1.0)
-    coefficients[2, :] =  (-1.0, -2.0,   1.0,  4.0,   1.0, -2.0, -1.0)
-    coefficients[3, :] =  (-1.0,  0.0,   3.0,  0.0,  -3.0,  0.0,  1.0)
-    coefficients[4, :] =  ( 1.0, -2.0,  -1.0,  4.0,  -1.0, -2.0,  1.0)
-    coefficients[5, :] =  (-1.0,  6.0, -15.0, 20.0, -15.0,  6.0, -1.0)
-
-    LAWSFilter= {'numKernels' : 6, 'kernelSize' : 7, 'filters' : 21,
-                 'coefficients': coefficients, 'names': names} 
-
-    return LAWSFilter
-
-def build_laws_masks(LAWSFilter):
-
-    """
-    masks = build_laws_masks(LAWSFilter)
-
-    takes the Laws Filter dictionary and builds the 21 7x7 kernel masks that are
-    used in Laws texture feature extraction. 
-
-    Parameters 
-    ----------
-
-    LAWSFilter : {dictionary}
-        a set of 6 length-7 Laws texture kernels
-
-    Returns 
-    ----------
-
-    masks : {list}
-        a list of 21 7x7 kernels (2D nd_array)
-
-    Examples:
-    use this along with FFT Pack to view the spatial frequency response. Create a 256x256 zero 
-    array and pad the first 7x7 with the Laws kernel and then get the 2D power spectrum and display
-    with pylab
-
-
-    LAWSFilter = build_laws_kernel()
-    mask = build_laws_masks(LAWSFilter)
-
-    mask_2 = masks[2]
-
-    z = NP.zeros(256*256, dtype=NP.float32).reshape(256, 256)
-    z[0:7, 0:7] = mask_0
-    x = abs(fftshift(fft2(z)))
-    pylab.imshow(x)
-
-    """
-
-    outer_indices = range(0, LAWSFilter['numKernels'])
-    mask_array = {}
-    count = 0
-    for i in outer_indices:
-        rowFilter = LAWSFilter['coefficients'][i]
-        colFilter = LAWSFilter['coefficients'][i]
-        matrix = NP.outer(rowFilter, colFilter)
-        mask_array[count] = 2.0*matrix
-        count = count + 1 
-        inner_indices = range(i+1, LAWSFilter['numKernels'])
-        for j in inner_indices:
-            colFilter = LAWSFilter['coefficients'][j]
-            matrix = NP.outer(rowFilter, colFilter) + NP.outer(colFilter, rowFilter)
-            mask_array[count] = matrix
-            count = count + 1 
-
-    return mask_array
-
-
-#
-#    test pattern generators for demo and test
-#
-
-def build_test_texture_discs():
-    """
-    discs = build_test_texture_discs()
-    
-    builds 4 discs with plane wave texture. used for test and demo
-
-    Parameters 
-    ----------
-
-    None
-
-    Returns 
-    ----------
-
-    discs : {nd_array}
-        a 512x512 image with 4 test discs (one per quadrant)
-
-    """
-    rows = 512
-    cols = 512
-    rad  = NP.pi / 180.0
-    test_image = NP.zeros(rows*cols, dtype=NP.float32).reshape(rows, cols)
-    [a, b] = NP.mgrid[0:rows, 0:cols]
-
-    test_image[0:255,0:255]     = NP.sin(4.0*rad*a[0:255,0:255])  + NP.sin(-4.0*rad*b[0:255,0:255]) 
-    test_image[256:511,256:511] = NP.sin(24.0*rad*a[0:255,0:255]) + NP.sin(20.0*rad*b[0:255,0:255])
-
-    test_image = test_image + test_image.min()
-    discs = build_test_unit_discs()
-    discs = discs * test_image 
-
-    return discs
-
-
-def build_test_discs():
-    """
-    test_image = build_test_discs()
-    build 4 discs of equal radius and different mean values for edge/blob testing
-    
-    Parameters 
-    ----------
-
-    None
-
-    Returns 
-    ----------
-
-    test_image : {nd_array}
-        a 512x512 image with 4 test discs (one per quadrant)
-
-    """
-    radius = 50
-    rows   = 512
-    cols   = 512
-    test_image = NP.zeros(rows*cols, dtype=NP.int16).reshape(rows, cols)
-    y_indices = NP.array(range(-radius, radius+1))
-    center_x = rows / 4
-    center_y = cols / 4
-
-    for i in y_indices:
-        x = math.sqrt(float(radius)**2 - float(i)**2)
-        # different raw mean levels
-        test_image[1*center_y+i, 1*center_x-x:1*center_x+x] = 80
-        test_image[1*center_y+i, 3*center_x-x:3*center_x+x] = 90
-        test_image[3*center_y+i, 1*center_x-x:1*center_x+x] = 100
-        test_image[3*center_y+i, 3*center_x-x:3*center_x+x] = 110
-
-    return test_image
-
-
-def build_test_unit_discs():
-    """
-    test_image = build_test_unit_discs()
-    build 2 discs of equal radius and same mean values for texture testing
-    
-    Parameters 
-    ----------
-
-    None
-
-    Returns 
-    ----------
-
-    test_image : {nd_array}
-        a 512x512 image with 4 test discs (one per quadrant)
-
-    """
-    radius = 50
-    rows   = 512
-    cols   = 512
-    test_image = NP.zeros(rows*cols, dtype=NP.int16).reshape(rows, cols)
-    y_indices = NP.array(range(-radius, radius+1))
-    center_x = rows / 4
-    center_y = cols / 4
-
-    for i in y_indices:
-        x = math.sqrt(float(radius)**2 - float(i)**2)
-        # different raw mean levels
-        test_image[1*center_y+i, 1*center_x-x:1*center_x+x] = 100
-        test_image[3*center_y+i, 3*center_x-x:3*center_x+x] = 100
-
-    return test_image
-
-
-def build_test_impulses():
-    """
-    test_image = build_test_impulses()
-
-    build 4 test impulses discs centered in the 4 discs. used
-    for testing filter kernels, esp. Laws' filter kernel. Filtering
-    with these test patterns will return Law's kernel outer product matrices.
-
-    Parameters 
-    ----------
-
-    None
-
-    Returns 
-    ----------
-
-    test_image : {nd_array}
-        a 512x512 image with 4 test discs (one per quadrant)
-
-    """
-    rows = 512
-    cols = 512
-    test_image = NP.zeros(rows*cols, dtype=NP.int16).reshape(rows, cols)
-    test_image[128,128] = 1 
-    test_image[378,128] = 1 
-    test_image[128,378] = 1 
-    test_image[378,378] = 1 
-
-    return test_image
-
-
-

Modified: trunk/scipy/ndimage/setup.py
===================================================================
--- trunk/scipy/ndimage/setup.py	2008-08-11 21:39:17 UTC (rev 4637)
+++ trunk/scipy/ndimage/setup.py	2008-08-11 23:08:56 UTC (rev 4638)
@@ -14,17 +14,6 @@
         include_dirs=['src']+[get_include()],
     )
 
-    config.add_extension('_segment',
-                         sources=['src/segment/Segmenter_EXT.c',
-                                  'src/segment/Segmenter_IMPL.c'],
-                         depends = ['src/segment/ndImage_Segmenter_structs.h']
-    )
-
-    config.add_extension('_register',
-                         sources=['src/register/Register_EXT.c',
-                                  'src/register/Register_IMPL.c']
-    )
-
     config.add_data_dir('tests')
 
     return config

Deleted: trunk/scipy/ndimage/tests/test_registration.py
===================================================================
--- trunk/scipy/ndimage/tests/test_registration.py	2008-08-11 21:39:17 UTC (rev 4637)
+++ trunk/scipy/ndimage/tests/test_registration.py	2008-08-11 23:08:56 UTC (rev 4638)
@@ -1,192 +0,0 @@
-import math
-import numpy as np
-import scipy.ndimage._registration as reg
-from numpy.testing import *
-
-def load_desc():
-    # this is for a 256x256x90 volume with 0.9375 x 0.9375 * 1.5 mm voxel sizes 
-    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 build_volume(imagedesc, S=[1500.0, 2500.0, 1000.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 = (255.0*np.exp(-(aa + bb + cc))).astype(np.uint8)
-
-    return volume3D
-
- # self.failUnless(diff(output, tcov) < eps)
-class TestRegistration(TestCase):
-
-    def test_affine_matrix_build_1(self):
-        "test_affine_matrix_build_1"
-        P = np.zeros(6)
-        M = reg.build_rotate_matrix(P)
-        E = np.eye(4)
-        match = (E==M).all()
-        assert_equal(match, True)
-        return
-
-    def test_affine_matrix_build_2(self):
-        "test_affine_matrix_build_2"
-        P = np.zeros(6)
-        P[0] = 1.0
-        M = reg.build_rotate_matrix(P)
-        E = np.array([
-                     [ 1. ,  0.        , 0.        , 0. ],
-                     [ 0. ,  0.9998477 , 0.01745241, 0. ],
-                     [ 0. , -0.01745241, 0.9998477 , 0. ],
-                     [ 0. ,  0.        , 0.        , 1. ]
-                     ])
-        assert_array_almost_equal(E, M, decimal=6)
-        return
-
-    def test_affine_matrix_build_3(self):
-        "test_affine_matrix_build_3"
-        P = np.zeros(6)
-        P[0] = 1.0
-        P[1] = 1.0
-        P[2] = 1.0
-        M = reg.build_rotate_matrix(P)
-        E = np.array([
-                     [ 0.99969541,  0.01744975,  0.01745241,  0. ],
-                     [-0.01775429,  0.9996901 ,  0.01744975,  0. ],
-                     [-0.0171425 , -0.01775429,  0.99969541,  0. ],
-                     [ 0.        ,  0.        ,  0.        ,  1. ]
-                     ])
-        assert_array_almost_equal(E, M, decimal=6)
-        return
-
-    @dec.slow
-    def test_autoalign_histogram_1(self):
-        "test_autoalign_histogram_1"
-        desc = load_desc()
-        gvol = build_volume(desc)
-        mat  = np.eye(4)
-        cost, joint = reg.check_alignment(gvol, mat, gvol, mat, ret_histo=1, lite=1)
-        # confirm that with lite=1 only have non-zero on the main diagonal
-        j_diag = joint.diagonal()
-        Z = np.diag(j_diag)
-        W = joint - Z
-        # clip the near-zero fuzz
-        W[abs(W) < 1e-10] = 0.0
-        assert_equal(W.max(), 0.0)
-        return
-
-    @dec.slow
-    def test_autoalign_histogram_2(self):
-        "test_autoalign_histogram_2"
-        desc = load_desc()
-        gvol = build_volume(desc)
-        mat  = np.eye(4)
-        cost, joint = reg.check_alignment(gvol, mat, gvol, mat, ret_histo=1, lite=0)
-        # confirm that with lite=0 DO have non-zero on the main diagonal
-        j_diag = joint.diagonal()
-        Z = np.diag(j_diag)
-        W = joint - Z
-        # clip the near-zero fuzz
-        W[abs(W) < 1e-10] = 0.0
-        s = (W.max() > 0.0)
-        # make sure there are off-diagonal values
-        assert_equal(s, True)
-        return
-
-    @dec.slow
-    def test_autoalign_ncc_value_1(self):
-        "test_autoalign_ncc_value_1"
-        desc = load_desc()
-        gvol = build_volume(desc)
-        mat  = np.eye(4)
-        cost = reg.check_alignment(gvol, mat, gvol, mat, method='ncc', lite=1)
-        # confirm the cross correlation is near 1.0 
-        t = abs(cost) + 0.0001
-        s = (t >= 1.0)
-        assert_equal(s, True)
-        return
-
-    @dec.slow
-    def test_autoalign_ncc_value_2(self):
-        "test_autoalign_ncc_value_2"
-        desc = load_desc()
-        gvol = build_volume(desc)
-        mat  = np.eye(4)
-        cost = reg.check_alignment(gvol, mat, gvol, mat, method='ncc', lite=0)
-        # confirm the cross correlation is near 1.0 
-        t = abs(cost) + 0.0001
-        s = (t >= 1.0)
-        assert_equal(s, True)
-        return
-
-    @dec.slow
-    def test_autoalign_nmi_value_1(self):
-        "test_autoalign_nmi_value_1"
-        desc = load_desc()
-        gvol = build_volume(desc)
-        mat  = np.eye(4)
-        cost = reg.check_alignment(gvol, mat, gvol, mat, method='nmi', lite=1)
-        # confirm the normalized mutual information is near -2.0 
-        assert_almost_equal(cost, -2.0, decimal=6)
-        return
-
-    @dec.slow
-    def test_autoalign_nmi_value_2(self):
-        "test_autoalign_nmi_value_2"
-        desc = load_desc()
-        gvol = build_volume(desc)
-        mat  = np.eye(4)
-        cost = reg.check_alignment(gvol, mat, gvol, mat, method='nmi', lite=0)
-        # confirm the normalized mutual information is near -2.0 
-        assert_almost_equal(cost, -1.7973048186515352, decimal=6)
-        return
-
-
-
-if __name__ == "__main__":
-    #nose.runmodule()
-    run_module_suite()
-
-
-

Deleted: trunk/scipy/ndimage/tests/test_segment.py
===================================================================
--- trunk/scipy/ndimage/tests/test_segment.py	2008-08-11 21:39:17 UTC (rev 4637)
+++ trunk/scipy/ndimage/tests/test_segment.py	2008-08-11 23:08:56 UTC (rev 4638)
@@ -1,159 +0,0 @@
-import math
-import numpy as NP
-import scipy.ndimage._segmenter as seg
-from numpy.testing import *
-
-def run_sobel():
-    img    = seg.build_test_discs()
-    filter = seg.build_2d_kernel(hiFilterCutoff=60.0)
-    fslice = seg.pre_filter(img, filter, low_threshold=0, high_threshold=255)
-    sobel_edge_image, sobel_stats = seg.sobel_image(fslice)
-    sobel_edge = seg.sobel_edges(sobel_edge_image, sobel_stats, sobel_threshold=0.8)
-    label_sobel, sobel_groups = seg.get_blobs(sobel_edge)
-    ROI      = seg.get_blob_regions(label_sobel, sobel_groups)
-    measures = seg.get_all_bounding_boxes(ROI)
-    thin_kernel = seg.build_morpho_thin_masks()
-    sobel_edges = seg.mat_filter(label_sobel, thin_kernel, ROI)
-    seg.get_voxel_measures(label_sobel, img, ROI)
-    means = ROI[:]['voxelMean']
-    return measures, means
-
-def run_canny():
-    img          = seg.build_test_discs()
-    filter       = seg.build_2d_kernel(hiFilterCutoff=60.0)
-    fslice       = seg.pre_filter(img, filter, low_threshold=0, high_threshold=255)
-    canny_kernel = seg.build_d_gauss_kernel()
-    horz, vert, imean = seg.canny_filter(fslice, canny_kernel)
-    mag, canny_stats  = seg.canny_nonmax_supress(horz, vert, imean)
-    canny_edge        = seg.canny_hysteresis(mag, canny_stats)
-    label_canny, canny_groups = seg.get_blobs(canny_edge)
-    ROI      = seg.get_blob_regions(label_canny, canny_groups)
-    measures = seg.get_all_bounding_boxes(ROI)
-    seg.get_voxel_measures(label_canny, img, ROI)
-    means = ROI[:]['voxelMean']
-    return measures, means
-
-def run_texture1():
-    filter      = seg.build_2d_kernel(hiFilterCutoff=60.0)
-    img         = seg.build_test_discs()
-    disc_mask   = seg.pre_filter(img, filter, low_threshold=50, high_threshold=255,
-                                 conv_binary=1)
-    label_disc_mask, disc_mask_groups = seg.get_blobs(disc_mask)
-    disc_ROI    = seg.get_blob_regions(label_disc_mask, disc_mask_groups)
-    laws_kernel = seg.build_laws_kernel() 
-    impulse     = seg.build_test_impulses()
-    calib       = seg.laws_texture_filter(impulse, label_disc_mask, laws_kernel,
-                                          ROI=disc_ROI, verbose=1)
-    kernels = calib[0]
-    x = laws_kernel['coefficients'][0]
-    m = NP.outer(x, x)
-    m = m * 2
-    laws_LL = kernels[0, 50-4:50+3, 50-3:50+4]
-    return m, laws_LL 
-
-
-def run_texture2():
-    filter = seg.build_2d_kernel(hiFilterCutoff=60.0)
-    img = seg.build_test_unit_discs()
-    disc = seg.pre_filter(img, filter, low_threshold=50, high_threshold=255)
-    disc_mask = seg.pre_filter(img, filter, low_threshold=50, high_threshold=255,
-                               conv_binary=1)
-    label_disc_mask, disc_mask_groups = seg.get_blobs(disc_mask)
-    disc_ROI = seg.get_blob_regions(label_disc_mask, disc_mask_groups)
-    laws_kernel = seg.build_laws_kernel() 
-    texture_img = seg.build_test_texture_discs()
-    seg.laws_texture_filter(texture_img, label_disc_mask, laws_kernel, ROI=disc_ROI,
-                            mean_feature=1, verbose=0)
-    tem = disc_ROI['TEM']
-    return tem 
-
-class TestSegment(TestCase):
-    def test_sobel(self):
-        # generate 4 discs, find the bounding boxes and 
-        # confirm the bounding boxes are at the true position
-        measures, voxel_means = run_sobel()
-        number = measures.size
-        _shortstruct = NP.dtype([('left', 'i'),
-                                 ('right', 'i'),
-                                 ('top', 'i'),
-                                 ('bottom', 'i')])
-
-        assert_equal(number, 4)
-        # load the ground truth
-        truth = NP.zeros(number, dtype=_shortstruct)
-        truth[0] = (76, 179, 179, 77) 
-        truth[1] = (332, 435, 179, 77)
-        truth[2] = (76, 179, 435, 333)
-        truth[3] = (332, 435, 435, 333)
-        match = (truth==measures).all()
-        assert_equal(match, True)
-        # load the ground truth for the bounding box test image mean value
-        voxel_truth = NP.zeros(number, dtype=NP.float64)
-        voxel_truth = (80.0, 90.0, 100.0, 110.0)
-        match = (voxel_truth==voxel_means).all()
-        assert_equal(match, True)
-
-        return
-
-    def test_canny(self):
-        # generate 4 discs, find the bounding boxes and 
-        # confirm the bounding boxes are at the true position
-        measures, voxel_means = run_canny()
-        number = measures.size
-        _shortstruct = NP.dtype([('left', 'i'),
-                                 ('right', 'i'),
-                                 ('top', 'i'),
-                                 ('bottom', 'i')])
-
-        assert_equal(number, 4)
-        # load the ground truth for the bounding box
-        truth = NP.zeros(number, dtype=_shortstruct)
-        truth[0] = (78,  177, 177, 79)
-        truth[1] = (334, 433, 177, 79)
-        truth[2] = (78,  177, 433, 335)
-        truth[3] = (334, 433, 433, 335)
-        match = (truth==measures).all()
-        assert_equal(match, True)
-        # load the ground truth for the bounding box test image mean value
-        voxel_truth = NP.zeros(number, dtype=NP.float64)
-        voxel_truth = (80.0, 90.0, 100.0, 110.0)
-        match = (voxel_truth==voxel_means).all()
-        assert_equal(match, True)
-
-        return
-
-    def test_texture1(self):
-        # [1] texture1 is delta functions and confirm the
-        #     filter result is outer product of the L kernel
-        M, Laws_LL = run_texture1()
-        match = (Laws_LL==M).all()
-        assert_equal(match, True)
-        return
-
-    def test_texture2(self):
-        # [2] texture2 is 2 plane waves and assert the 20-element feature
-        #     vector for each disc is correct
-        tem = run_texture2()
-        tem0 = tem[0]
-        tem1 = tem[1]
-        truth_tem0 = NP.array(
-                        [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
-                          0.        ,  0.13306101,  0.08511007,  0.05084148,  0.07550675,
-                          0.4334695 ,  0.03715914,  0.00289055,  0.02755581,  0.48142046,
-                          0.03137803,  0.00671277,  0.51568902,  0.01795249,  0.49102375,  1.
-                        ], dtype=NP.float32)
-        truth_tem1 = NP.array(
-                        [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
-                          0.        ,  0.02970393,  0.00164266,  0.00922416,  0.01221788,
-                          0.51485199,  0.03298925,  0.02212243,  0.01912871,  0.48350537,
-                          0.01125561,  0.00826189,  0.49437219,  0.00526817,  0.49736592,  1.
-                        ], dtype=NP.float32)
-
-        assert_array_almost_equal(tem0, truth_tem0, decimal=6)
-        assert_array_almost_equal(tem1, truth_tem1, decimal=6)
-
-        return
-
-if __name__ == "__main__":
-    nose.runmodule()
-



More information about the Scipy-svn mailing list