[Scipy-svn] r6339 - in trunk/scipy/ndimage: . src tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sat Apr 24 16:59:32 EDT 2010
Author: stefan
Date: 2010-04-24 15:59:32 -0500 (Sat, 24 Apr 2010)
New Revision: 6339
Modified:
trunk/scipy/ndimage/measurements.py
trunk/scipy/ndimage/setup.py
trunk/scipy/ndimage/src/nd_image.c
trunk/scipy/ndimage/tests/test_ndimage.py
Log:
Merge branch 'ndimage_measurements_rewrite' of http://broad.mit.edu/~thouis/scipy into measurements
Modified: trunk/scipy/ndimage/measurements.py
===================================================================
--- trunk/scipy/ndimage/measurements.py 2010-04-23 00:21:57 UTC (rev 6338)
+++ trunk/scipy/ndimage/measurements.py 2010-04-24 20:59:32 UTC (rev 6339)
@@ -34,6 +34,7 @@
import _ni_support
import _nd_image
import morphology
+import time
def label(input, structure = None, output = None):
"""
@@ -180,7 +181,160 @@
max_label = input.max()
return _nd_image.find_objects(input, max_label)
-def sum(input, labels=None, index=None):
+def labeled_comprehension(input, labels, index, func, out_dtype, default, pass_positions=False):
+ '''Roughly equivalent to [func(input[labels == i]) for i in index].
+
+ Special cases:
+ - index a scalar: returns a single value
+ - index is None: returns func(inputs[labels > 0])
+
+ func will be called with linear indices as a second argument if
+ pass_positions is True.
+ '''
+
+ as_scalar = numpy.isscalar(index)
+ input = numpy.asarray(input)
+
+ if pass_positions:
+ positions = numpy.arange(input.size).reshape(input.shape)
+
+ if labels is None:
+ if index is not None:
+ raise ValueError, "index without defined labels"
+ if not pass_positions:
+ return func(input.ravel())
+ else:
+ return func(input.ravel(), positions.ravel())
+
+ try:
+ input, labels = numpy.broadcast_arrays(input, labels)
+ except ValueError:
+ raise ValueError, "input and labels must have the same shape (excepting dimensions with width 1)"
+
+ if index is None:
+ if not pass_positions:
+ return func(input[labels > 0])
+ else:
+ return func(input[labels > 0], positions[labels > 0])
+
+ index = numpy.atleast_1d(index)
+ if any(index.astype(labels.dtype).astype(index.dtype) != index):
+ raise ValueError, "Cannot convert index values from <%s> to <%s> (labels' type) without loss of precision"%(index.dtype, labels.dtype)
+ index = index.astype(labels.dtype)
+
+ # optimization: find min/max in index, and select those parts of labels, input, and positions
+ lo = index.min()
+ hi = index.max()
+ mask = (labels >= lo) & (labels <= hi)
+
+ # this also ravels the arrays
+ labels = labels[mask]
+ input = input[mask]
+ if pass_positions:
+ positions = positions[mask]
+
+ # sort everything by labels
+ label_order = labels.argsort()
+ labels = labels[label_order]
+ input = input[label_order]
+ if pass_positions:
+ positions = positions[label_order]
+
+ index_order = index.argsort()
+ sorted_index = index[index_order]
+
+ def do_map(inputs, output):
+ '''labels must be sorted'''
+
+ nlabels = labels.size
+ nidx = sorted_index.size
+
+ # Find boundaries for each stretch of constant labels
+ # This could be faster, but we already paid N log N to sort labels.
+ lo = numpy.searchsorted(labels, sorted_index, side='left')
+ hi = numpy.searchsorted(labels, sorted_index, side='right')
+
+ for i, l, h in zip(range(nidx), lo, hi):
+ if l == h:
+ continue
+ idx = sorted_index[i]
+ output[i] = func(*[inp[l:h] for inp in inputs])
+
+ temp = numpy.empty(index.shape, out_dtype)
+ temp[:] = default
+ if not pass_positions:
+ do_map([input], temp)
+ else:
+ do_map([input, positions], temp)
+ output = numpy.zeros(index.shape, out_dtype)
+ output[index_order] = temp
+
+ if as_scalar:
+ output = output[0]
+
+ return output
+
+def _stats(input, labels = None, index = None, do_sum2=False):
+ '''returns count, sum, and optionally sum^2 by label'''
+
+ def single_group(vals):
+ if do_sum2:
+ return vals.size, vals.sum(), (vals * vals.conjugate()).sum()
+ else:
+ return vals.size, vals.sum()
+
+ if labels is None:
+ return single_group(input)
+
+ # ensure input and labels match sizes
+ input, labels = numpy.broadcast_arrays(input, labels)
+
+ if index is None:
+ return single_group(input[labels > 0])
+
+ if numpy.isscalar(index):
+ return single_group(input[labels == index])
+
+ # remap labels to unique integers if necessary, or if the largest
+ # label is larger than the number of values.
+ if ((not numpy.issubdtype(labels.dtype, numpy.int)) or
+ (labels.min() < 0) or (labels.max() > labels.size)):
+ unique_labels, new_labels = numpy.unique1d(labels, return_inverse=True)
+
+ counts = numpy.bincount(new_labels)
+ sums = numpy.bincount(new_labels, weights=input.ravel())
+ if do_sum2:
+ sums2 = numpy.bincount(new_labels, weights=(input * input.conjugate()).ravel())
+
+ idxs = numpy.searchsorted(unique_labels, index)
+ # make all of idxs valid
+ idxs[idxs >= unique_labels.size] = 0
+ found = (unique_labels[idxs] == index)
+ else:
+ # labels are an integer type, and there aren't too many, so
+ # call bincount directly.
+ counts = numpy.bincount(labels.ravel())
+ sums = numpy.bincount(labels.ravel(), weights=input.ravel())
+ if do_sum2:
+ sums2 = numpy.bincount(labels.ravel(), weights=(input * input.conjugate()).ravel())
+
+ # make sure all index values are valid
+ idxs = numpy.asanyarray(index, numpy.int).copy()
+ found = (idxs >= 0) & (idxs < counts.size)
+ idxs[~ found] = 0
+
+ counts = counts[idxs]
+ counts[~ found] = 0
+ sums = sums[idxs]
+ sums[~ found] = 0
+ if not do_sum2:
+ return (counts, sums)
+ sums2 = sums2[idxs]
+ sums2[~ found] = 0
+ return (counts, sums, sums2)
+
+
+def sum(input, labels = None, index = None):
"""Calculate the sum of the values of the array.
:Parameters:
@@ -201,248 +355,253 @@
[1.0, 5.0]
"""
- input = numpy.asarray(input)
- if numpy.iscomplexobj(input):
- raise TypeError, 'Complex type not supported'
- if labels is not None:
- labels = numpy.asarray(labels)
- labels = _broadcast(labels, input.shape)
+ count, sum = _stats(input, labels, index)
+ return sum
- if labels.shape != input.shape:
- raise RuntimeError, 'input and labels shape are not equal'
- if index is not None:
- T = getattr(index,'dtype',numpy.int32)
- if T not in [numpy.int8, numpy.int16, numpy.int32,
- numpy.uint8, numpy.uint16, numpy.bool]:
- raise ValueError("Invalid index type")
- index = numpy.asarray(index,dtype=T)
- return _nd_image.statistics(input, labels, index, 0)
+def mean(input, labels = None, index = None):
+ """Calculate the mean of the values of an array at labels.
+ Labels must be None or an array that can be broadcast to the input.
-def mean(input, labels = None, index = None):
- """Calculate the mean of the values of the array.
-
- The index parameter is a single label number or a sequence of
- label numbers of the objects to be measured. If index is None, all
- values are used where labels is larger than zero.
+ Index must be None, a single label or sequence of labels. If
+ None, the mean for all values where label is greater than 0 is
+ calculated.
"""
- input = numpy.asarray(input)
- if numpy.iscomplexobj(input):
- raise TypeError, 'Complex type not supported'
- if labels is not None:
- labels = numpy.asarray(labels)
- labels = _broadcast(labels, input.shape)
- if labels.shape != input.shape:
- raise RuntimeError, 'input and labels shape are not equal'
- return _nd_image.statistics(input, labels, index, 1)
+ count, sum = _stats(input, labels, index)
+ return sum / numpy.asanyarray(count).astype(numpy.float)
-
def variance(input, labels = None, index = None):
- """Calculate the variance of the values of the array.
+ """Calculate the variance of the values of an array at labels.
- The index parameter is a single label number or a sequence of
- label numbers of the objects to be measured. If index is None, all
- values are used where labels is larger than zero.
+ Labels must be None or an array of the same dimensions as the input.
+
+ Index must be None, a single label or sequence of labels. If
+ none, all values where label is greater than zero are used.
"""
- input = numpy.asarray(input)
- if numpy.iscomplexobj(input):
- raise TypeError, 'Complex type not supported'
- if labels is not None:
- labels = numpy.asarray(labels)
- labels = _broadcast(labels, input.shape)
- if labels.shape != input.shape:
- raise RuntimeError, 'input and labels shape are not equal'
- return _nd_image.statistics(input, labels, index, 2)
+ count, sum, sum2 = _stats(input, labels, index, do_sum2=True)
+ mean = sum / numpy.asanyarray(count).astype(numpy.float)
+ mean2 = sum2 / numpy.asanyarray(count).astype(numpy.float)
+ return mean2 - (mean * mean.conjugate())
def standard_deviation(input, labels = None, index = None):
- """Calculate the standard deviation of the values of the array.
+ """Calculate the standard deviation of the values of an array at labels.
- The index parameter is a single label number or a sequence of
- label numbers of the objects to be measured. If index is None, all
- values are used where labels is larger than zero.
+ Labels must be None or an array of the same dimensions as the input.
+
+ Index must be None, a single label or sequence of labels. If
+ none, all values where label is greater than zero are used.
"""
- var = variance(input, labels, index)
- if (isinstance(var, types.ListType)):
- return [math.sqrt(x) for x in var]
- else:
- return math.sqrt(var)
+ return numpy.sqrt(variance(input, labels, index))
-def minimum(input, labels = None, index = None):
- """Calculate the minimum of the values of the array.
+def _select(input, labels = None, index = None, find_min=False, find_max=False, find_min_positions=False, find_max_positions=False):
+ '''returns min, max, or both, plus positions if requested'''
- The index parameter is a single label number or a sequence of
- label numbers of the objects to be measured. If index is None, all
- values are used where labels is larger than zero.
- """
- input = numpy.asarray(input)
- if numpy.iscomplexobj(input):
- raise TypeError, 'Complex type not supported'
- if labels is not None:
- labels = numpy.asarray(labels)
- labels = _broadcast(labels, input.shape)
- if labels.shape != input.shape:
- raise RuntimeError, 'input and labels shape are not equal'
- return _nd_image.statistics(input, labels, index, 3)
+ find_positions = find_min_positions or find_max_positions
+ positions = None
+ if find_positions:
+ positions = numpy.arange(input.size).reshape(input.shape)
+ def single_group(vals, positions):
+ result = []
+ if find_min:
+ result += [vals.min()]
+ if find_min_positions:
+ result += [positions[vals == vals.min()][0]]
+ if find_max:
+ result += [vals.max()]
+ if find_max_positions:
+ result += [positions[vals == vals.max()][0]]
+ return result
+
+ if labels is None:
+ return single_group(input, positions)
-def maximum(input, labels=None, index=None):
- """Return the maximum input value.
+ # ensure input and labels match sizes
+ input, labels = numpy.broadcast_arrays(input, labels)
- The index parameter is a single label number or a sequence of
- label numbers of the objects to be measured. If index is None, all
- values are used where labels is larger than zero.
+ if index is None:
+ mask = (labels > 0)
+ return single_group(input[mask], positions[mask] if find_positions else None)
- """
- input = numpy.asarray(input)
- if numpy.iscomplexobj(input):
- raise TypeError, 'Complex type not supported'
- if labels is not None:
- labels = numpy.asarray(labels)
- labels = _broadcast(labels, input.shape)
+ if numpy.isscalar(index):
+ mask = (labels == index)
+ return single_group(input[mask], positions[mask] if find_positions else None)
- if labels.shape != input.shape:
- raise RuntimeError, 'input and labels shape are not equal'
- return _nd_image.statistics(input, labels, index, 4)
+ order = input.ravel().argsort()
+ input = input.ravel()[order]
+ labels = labels.ravel()[order]
+ if find_positions:
+ positions = positions.ravel()[order]
+ # remap labels to unique integers if necessary, or if the largest
+ # label is larger than the number of values.
+ if ((not numpy.issubdtype(labels.dtype, numpy.int)) or
+ (labels.min() < 0) or (labels.max() > labels.size)):
+ # remap labels, and indexes
+ unique_labels, labels = numpy.unique1d(labels, return_inverse=True)
+ idxs = numpy.searchsorted(unique_labels, index)
-def _index_to_position(index, shape):
- """Convert a linear index to a position"""
- if len(shape) > 0:
- pos = []
- stride = numpy.multiply.reduce(shape)
- for size in shape:
- stride = stride // size
- pos.append(index // stride)
- index -= pos[-1] * stride
- return tuple(pos)
+ # make all of idxs valid
+ idxs[idxs >= unique_labels.size] = 0
+ found = (unique_labels[idxs] == index)
else:
- return 0
+ # labels are an integer type, and there aren't too many.
+ idxs = numpy.asanyarray(index, numpy.int).copy()
+ found = (idxs >= 0) & (idxs <= labels.max())
+
+ idxs[~ found] = labels.max() + 1
+ result = []
+ if find_min:
+ mins = numpy.zeros(labels.max() + 2, input.dtype)
+ mins[labels[::-1]] = input[::-1]
+ result += [mins[idxs]]
+ if find_min_positions:
+ minpos = numpy.zeros(labels.max() + 2)
+ minpos[labels[::-1]] = positions[::-1]
+ result += [minpos[idxs]]
+ if find_max:
+ maxs = numpy.zeros(labels.max() + 2, input.dtype)
+ maxs[labels] = input
+ result += [maxs[idxs]]
+ if find_max_positions:
+ maxpos = numpy.zeros(labels.max() + 2)
+ maxpos[labels] = positions
+ result += [maxpos[idxs]]
+ return result
+def minimum(input, labels = None, index = None):
+ """Calculate the minimum of the values of an array at labels.
+
+ Labels must be None or an array of the same dimensions as the input.
+
+ Index must be None, a single label or sequence of labels. If
+ none, all values where label is greater than zero are used.
+ """
+ return _select(input, labels, index, find_min=True)[0]
+
+def maximum(input, labels = None, index = None):
+ """Calculate the maximum of the values of an array at labels.
+
+ Labels must be None or an array of the same dimensions as the input.
+
+ Index must be None, a single label or sequence of labels. If
+ none, all values where label is greater than zero are used.
+ """
+ return _select(input, labels, index, find_max=True)[0]
+
def minimum_position(input, labels = None, index = None):
- """Find the position of the minimum of the values of the array.
+ """Find the positions of the minimums of the values of an array at labels.
- The index parameter is a single label number or a sequence of
- label numbers of the objects to be measured. If index is None, all
- values are used where labels is larger than zero.
+ Labels must be None or an array of the same dimensions as the input.
+
+ Index must be None, a single label or sequence of labels. If
+ none, all values where label is greater than zero are used.
"""
- input = numpy.asarray(input)
- if numpy.iscomplexobj(input):
- raise TypeError, 'Complex type not supported'
- if labels is not None:
- labels = numpy.asarray(labels)
- labels = _broadcast(labels, input.shape)
+
+ dims = numpy.array(numpy.asarray(input).shape)
+ # see numpy.unravel_index to understand this line.
+ dim_prod = numpy.cumprod([1] + list(dims[:0:-1]))[::-1]
- if labels.shape != input.shape:
- raise RuntimeError, 'input and labels shape are not equal'
- pos = _nd_image.statistics(input, labels, index, 5)
- if (isinstance(pos, types.ListType)):
- return [_index_to_position(x, input.shape) for x in pos]
- else:
- return _index_to_position(pos, input.shape)
+ result = _select(input, labels, index, find_min_positions=True)[0]
+ if numpy.isscalar(result):
+ return tuple((result // dim_prod) % dims)
+ return [tuple(v) for v in (result.reshape(-1, 1) // dim_prod) % dims]
+
def maximum_position(input, labels = None, index = None):
- """Find the position of the maximum of the values of the array.
+ """Find the positions of the maximums of the values of an array at labels.
- The index parameter is a single label number or a sequence of
- label numbers of the objects to be measured. If index is None, all
- values are used where labels is larger than zero.
+ Labels must be None or an array of the same dimensions as the input.
+
+ Index must be None, a single label or sequence of labels. If
+ none, all values where label is greater than zero are used.
"""
- input = numpy.asarray(input)
- if numpy.iscomplexobj(input):
- raise TypeError, 'Complex type not supported'
- if labels is not None:
- labels = numpy.asarray(labels)
- labels = _broadcast(labels, input.shape)
+
+ dims = numpy.array(numpy.asarray(input).shape)
+ # see numpy.unravel_index to understand this line.
+ dim_prod = numpy.cumprod([1] + list(dims[:0:-1]))[::-1]
- if labels.shape != input.shape:
- raise RuntimeError, 'input and labels shape are not equal'
- pos = _nd_image.statistics(input, labels, index, 6)
- if (isinstance(pos, types.ListType)):
- return [_index_to_position(x, input.shape) for x in pos]
- else:
- return _index_to_position(pos, input.shape)
+ result = _select(input, labels, index, find_max_positions=True)[0]
+ if numpy.isscalar(result):
+ return tuple((result // dim_prod) % dims)
+ return [tuple(v) for v in (result.reshape(-1, 1) // dim_prod) % dims]
+
def extrema(input, labels = None, index = None):
- """Calculate the minimum, the maximum and their positions of the
- values of the array.
+ """Calculate the minimums and maximums of the values of an array
+ at labels, along with their positions.
- The index parameter is a single label number or a sequence of
- label numbers of the objects to be measured. If index is None, all
- values are used where labels is larger than zero.
+ Labels must be None or an array of the same dimensions as the input.
+
+ Index must be None, a single label or sequence of labels. If
+ none, all values where label is greater than zero are used.
+
+ Returns: minimums, maximums, min_positions, max_positions
"""
- input = numpy.asarray(input)
- if numpy.iscomplexobj(input):
- raise TypeError, 'Complex type not supported'
- if labels is not None:
- labels = numpy.asarray(labels)
- labels = _broadcast(labels, input.shape)
+
+ dims = numpy.array(numpy.asarray(input).shape)
+ # see numpy.unravel_index to understand this line.
+ dim_prod = numpy.cumprod([1] + list(dims[:0:-1]))[::-1]
- if labels.shape != input.shape:
- raise RuntimeError, 'input and labels shape are not equal'
+ minimums, min_positions, maximums, max_positions = _select(input, labels, index,
+ find_min=True, find_max=True,
+ find_min_positions=True, find_max_positions=True)
- min, max, minp, maxp = _nd_image.statistics(input, labels, index, 7)
- if (isinstance(minp, types.ListType)):
- minp = [_index_to_position(x, input.shape) for x in minp]
- maxp = [_index_to_position(x, input.shape) for x in maxp]
- else:
- minp = _index_to_position(minp, input.shape)
- maxp = _index_to_position(maxp, input.shape)
- return min, max, minp, maxp
+ if numpy.isscalar(minimums):
+ return minimums, maximums, tuple((min_positions // dim_prod) % dims), tuple((max_positions // dim_prod) % dims)
+ min_positions = [tuple(v) for v in (min_positions.reshape(-1, 1) // dim_prod) % dims]
+ max_positions = [tuple(v) for v in (max_positions.reshape(-1, 1) // dim_prod) % dims]
+ return minimums, maximums, min_positions, max_positions
+
def center_of_mass(input, labels = None, index = None):
- """Calculate the center of mass of of the array.
+ """Calculate the center of mass of the values of an array at labels.
- The index parameter is a single label number or a sequence of
- label numbers of the objects to be measured. If index is None, all
- values are used where labels is larger than zero.
+ Labels must be None or an array of the same dimensions as the input.
+
+ Index must be None, a single label or sequence of labels. If
+ none, all values where label is greater than zero are used.
"""
- input = numpy.asarray(input)
- if numpy.iscomplexobj(input):
- raise TypeError, 'Complex type not supported'
- if labels is not None:
- labels = numpy.asarray(labels)
- labels = _broadcast(labels, input.shape)
- if labels.shape != input.shape:
- raise RuntimeError, 'input and labels shape are not equal'
- return _nd_image.center_of_mass(input, labels, index)
+ normalizer = sum(input, labels, index)
+ grids = numpy.ogrid[[slice(0, i) for i in input.shape]]
+ results = [sum(input * grids[dir].astype(float), labels, index) / normalizer for dir in range(input.ndim)]
+
+ if numpy.isscalar(results[0]):
+ return tuple(results)
+ return [tuple(v) for v in numpy.array(results).T]
+
def histogram(input, min, max, bins, labels = None, index = None):
- """Calculate a histogram of of the array.
+ """Calculate the histogram of the values of an array at labels.
- The histogram is defined by its minimum and maximum value and the
+ Labels must be None or an array of the same dimensions as the input.
+
+ The histograms are defined by the minimum and maximum values and the
number of bins.
- The index parameter is a single label number or a sequence of
- label numbers of the objects to be measured. If index is None, all
- values are used where labels is larger than zero.
+ Index must be None, a single label or sequence of labels. If
+ none, all values where label is greater than zero are used.
"""
- input = numpy.asarray(input)
- if numpy.iscomplexobj(input):
- raise TypeError, 'Complex type not supported'
- if labels is not None:
- labels = numpy.asarray(labels)
- labels = _broadcast(labels, input.shape)
+
+ _bins = numpy.linspace(min, max, bins + 1)
- if labels.shape != input.shape:
- raise RuntimeError, 'input and labels shape are not equal'
- if bins < 1:
- raise RuntimeError, 'number of bins must be >= 1'
- if min >= max:
- raise RuntimeError, 'min must be < max'
- return _nd_image.histogram(input, min, max, bins, labels, index)
+ def _hist(vals):
+ return numpy.histogram(vals, _bins)[0]
+ return labeled_comprehension(input, labels, index, _hist, object, None, pass_positions=False)
+
def watershed_ift(input, markers, structure = None, output = None):
"""Apply watershed from markers using a iterative forest transform
algorithm.
@@ -489,23 +648,3 @@
output, return_value = _ni_support._get_output(output, input)
_nd_image.watershed_ift(input, markers, structure, output)
return return_value
-
-def _broadcast(arr, sshape):
- """Return broadcast view of arr, else return None."""
- ashape = arr.shape
- return_value = numpy.zeros(sshape, arr.dtype)
- # Just return arr if they have the same shape
- if sshape == ashape:
- return arr
- srank = len(sshape)
- arank = len(ashape)
-
- aslices = []
- sslices = []
- for i in range(arank):
- aslices.append(slice(0, ashape[i], 1))
-
- for i in range(srank):
- sslices.append(slice(0, sshape[i], 1))
- return_value[sslices] = arr[aslices]
- return return_value
Modified: trunk/scipy/ndimage/setup.py
===================================================================
--- trunk/scipy/ndimage/setup.py 2010-04-23 00:21:57 UTC (rev 6338)
+++ trunk/scipy/ndimage/setup.py 2010-04-24 20:59:32 UTC (rev 6339)
@@ -12,6 +12,7 @@
"src/ni_measure.c",
"src/ni_morphology.c","src/ni_support.c"],
include_dirs=['src']+[get_include()],
+ extra_compile_args=['-Wall'],
)
config.add_data_dir('tests')
Modified: trunk/scipy/ndimage/src/nd_image.c
===================================================================
--- trunk/scipy/ndimage/src/nd_image.c 2010-04-23 00:21:57 UTC (rev 6338)
+++ trunk/scipy/ndimage/src/nd_image.c 2010-04-24 20:59:32 UTC (rev 6339)
@@ -731,399 +731,6 @@
return PyErr_Occurred() ? NULL : Py_BuildValue("");
}
-static int _NI_GetIndices(PyObject* indices_object,
- maybelong** result_indices, maybelong* min_label,
- maybelong* max_label, maybelong* n_results)
-{
- maybelong *indices = NULL, n_indices, ii;
-
- if (indices_object == Py_None) {
- *min_label = -1;
- *n_results = 1;
- } else {
- n_indices = NI_ObjectToLongSequenceAndLength(indices_object, &indices);
- if (n_indices < 0)
- goto exit;
- if (n_indices < 1) {
- PyErr_SetString(PyExc_RuntimeError, "no correct indices provided");
- goto exit;
- } else {
- *min_label = *max_label = indices[0];
- if (*min_label < 0) {
- PyErr_SetString(PyExc_RuntimeError,
- "negative indices not allowed");
- goto exit;
- }
- for(ii = 1; ii < n_indices; ii++) {
- if (indices[ii] < 0) {
- PyErr_SetString(PyExc_RuntimeError,
- "negative indices not allowed");
- goto exit;
- }
- if (indices[ii] < *min_label)
- *min_label = indices[ii];
- if (indices[ii] > *max_label)
- *max_label = indices[ii];
- }
- *result_indices = (maybelong*)malloc((*max_label - *min_label + 1) *
- sizeof(maybelong));
- if (!*result_indices) {
- PyErr_NoMemory();
- goto exit;
- }
- for(ii = 0; ii < *max_label - *min_label + 1; ii++)
- (*result_indices)[ii] = -1;
- *n_results = 0;
- for(ii = 0; ii < n_indices; ii++) {
- if ((*result_indices)[indices[ii] - *min_label] >= 0) {
- PyErr_SetString(PyExc_RuntimeError, "duplicate index");
- goto exit;
- }
- (*result_indices)[indices[ii] - *min_label] = ii;
- ++(*n_results);
- }
- }
- }
- exit:
- if (indices)
- free(indices);
- return PyErr_Occurred() == NULL;
-}
-
-
-PyObject* _NI_BuildMeasurementResultArrayObject(maybelong n_results,
- PyArrayObject** values)
-{
- PyObject *result = NULL;
- if (n_results > 1) {
- result = PyList_New(n_results);
- if (result) {
- maybelong ii;
- for(ii = 0; ii < n_results; ii++) {
- PyList_SET_ITEM(result, ii, (PyObject*)values[ii]);
- Py_XINCREF(values[ii]);
- }
- }
- } else {
- result = (PyObject*)values[0];
- Py_XINCREF(values[0]);
- }
- return result;
-}
-
-
-PyObject* _NI_BuildMeasurementResultDouble(maybelong n_results,
- double* values)
-{
- PyObject *result = NULL;
- if (n_results > 1) {
- result = PyList_New(n_results);
- if (result) {
- int ii;
- for(ii = 0; ii < n_results; ii++) {
- PyObject* val = PyFloat_FromDouble(values[ii]);
- if (!val) {
- Py_XDECREF(result);
- return NULL;
- }
- PyList_SET_ITEM(result, ii, val);
- }
- }
- } else {
- result = Py_BuildValue("d", values[0]);
- }
- return result;
-}
-
-
-PyObject* _NI_BuildMeasurementResultDoubleTuple(maybelong n_results,
- int tuple_size, double* values)
-{
- PyObject *result = NULL;
- maybelong ii;
- int jj;
-
- if (n_results > 1) {
- result = PyList_New(n_results);
- if (result) {
- for(ii = 0; ii < n_results; ii++) {
- PyObject* val = PyTuple_New(tuple_size);
- if (!val) {
- Py_XDECREF(result);
- return NULL;
- }
- for(jj = 0; jj < tuple_size; jj++) {
- maybelong idx = jj + ii * tuple_size;
- PyTuple_SetItem(val, jj, PyFloat_FromDouble(values[idx]));
- if (PyErr_Occurred()) {
- Py_XDECREF(result);
- return NULL;
- }
- }
- PyList_SET_ITEM(result, ii, val);
- }
- }
- } else {
- result = PyTuple_New(tuple_size);
- if (result) {
- for(ii = 0; ii < tuple_size; ii++) {
- PyTuple_SetItem(result, ii, PyFloat_FromDouble(values[ii]));
- if (PyErr_Occurred()) {
- Py_XDECREF(result);
- return NULL;
- }
- }
- }
- }
- return result;
-}
-
-
-PyObject* _NI_BuildMeasurementResultInt(maybelong n_results,
- maybelong* values)
-{
- PyObject *result = NULL;
- if (n_results > 1) {
- result = PyList_New(n_results);
- if (result) {
- maybelong ii;
- for(ii = 0; ii < n_results; ii++) {
- PyObject* val = PyInt_FromLong(values[ii]);
- if (!val) {
- Py_XDECREF(result);
- return NULL;
- }
- PyList_SET_ITEM(result, ii, val);
- }
- }
- } else {
- result = Py_BuildValue("l", values[0]);
- }
- return result;
-}
-
-
-static PyObject *Py_Statistics(PyObject *obj, PyObject *args)
-{
- PyArrayObject *input = NULL, *labels = NULL;
- PyObject *indices_object, *result = NULL;
- PyObject *res1 = NULL, *res2 = NULL, *res3 = NULL, *res4 = NULL;
- double *dresult1 = NULL, *dresult2 = NULL;
- maybelong *lresult1 = NULL, *lresult2 = NULL;
- maybelong min_label, max_label, *result_indices = NULL, n_results, ii;
- int type;
-
- if (!PyArg_ParseTuple(args, "O&O&Oi", NI_ObjectToInputArray, &input,
- NI_ObjectToOptionalInputArray, &labels, &indices_object, &type))
- goto exit;
-
- if (!_NI_GetIndices(indices_object, &result_indices, &min_label,
- &max_label, &n_results))
- goto exit;
-
- if (type >= 0 && type <= 7) {
- dresult1 = (double*)malloc(n_results * sizeof(double));
- if (!dresult1) {
- PyErr_NoMemory();
- goto exit;
- }
- }
- if (type == 2 || type == 7) {
- dresult2 = (double*)malloc(n_results * sizeof(double));
- if (!dresult2) {
- PyErr_NoMemory();
- goto exit;
- }
- }
- if (type == 1 || type == 2 || (type >= 5 && type <= 7)) {
- lresult1 = (maybelong*)malloc(n_results * sizeof(maybelong));
- if (!lresult1) {
- PyErr_NoMemory();
- goto exit;
- }
- }
- if (type == 7) {
- lresult2 = (maybelong*)malloc(n_results * sizeof(maybelong));
- if (!lresult2) {
- PyErr_NoMemory();
- goto exit;
- }
- }
- switch(type) {
- case 0:
- if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
- n_results, dresult1, NULL, NULL, NULL, NULL, NULL, NULL))
- goto exit;
- result = _NI_BuildMeasurementResultDouble(n_results, dresult1);
- break;
- case 1:
- if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
- n_results, dresult1, lresult1, NULL, NULL, NULL, NULL, NULL))
- goto exit;
- for(ii = 0; ii < n_results; ii++)
- dresult1[ii] = lresult1[ii] > 0 ? dresult1[ii] / lresult1[ii] : 0.0;
-
- result = _NI_BuildMeasurementResultDouble(n_results, dresult1);
- break;
- case 2:
- if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
- n_results, dresult1, lresult1, dresult2, NULL, NULL, NULL, NULL))
- goto exit;
- result = _NI_BuildMeasurementResultDouble(n_results, dresult2);
- break;
- case 3:
- if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
- n_results, NULL, NULL, NULL, dresult1, NULL, NULL, NULL))
- goto exit;
- result = _NI_BuildMeasurementResultDouble(n_results, dresult1);
- break;
- case 4:
- if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
- n_results, NULL, NULL, NULL, NULL, dresult1, NULL, NULL))
- goto exit;
- result = _NI_BuildMeasurementResultDouble(n_results, dresult1);
- break;
- case 5:
- if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
- n_results, NULL, NULL, NULL, dresult1, NULL, lresult1, NULL))
- goto exit;
- result = _NI_BuildMeasurementResultInt(n_results, lresult1);
- break;
- case 6:
- if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
- n_results, NULL, NULL, NULL, NULL, dresult1, NULL, lresult1))
- goto exit;
- result = _NI_BuildMeasurementResultInt(n_results, lresult1);
- break;
- case 7:
- if (!NI_Statistics(input, labels, min_label, max_label, result_indices,
- n_results, NULL, NULL, NULL, dresult1, dresult2,
- lresult1, lresult2))
- goto exit;
- res1 = _NI_BuildMeasurementResultDouble(n_results, dresult1);
- res2 = _NI_BuildMeasurementResultDouble(n_results, dresult2);
- res3 = _NI_BuildMeasurementResultInt(n_results, lresult1);
- res4 = _NI_BuildMeasurementResultInt(n_results, lresult2);
- if (!res1 || !res2 || !res3 || !res4)
- goto exit;
- result = Py_BuildValue("OOOO", res1, res2, res3, res4);
- break;
- default:
- PyErr_SetString(PyExc_RuntimeError, "operation not supported");
- goto exit;
- }
-
- exit:
- Py_XDECREF(input);
- Py_XDECREF(labels);
- if (result_indices)
- free(result_indices);
- if (dresult1)
- free(dresult1);
- if (dresult2)
- free(dresult2);
- if (lresult1)
- free(lresult1);
- if (lresult2)
- free(lresult2);
- return result;
-}
-
-
-static PyObject *Py_CenterOfMass(PyObject *obj, PyObject *args)
-{
- PyArrayObject *input = NULL, *labels = NULL;
- PyObject *indices_object, *result = NULL;
- double *center_of_mass = NULL;
- maybelong min_label, max_label, *result_indices = NULL, n_results;
-
- if (!PyArg_ParseTuple(args, "O&O&O", NI_ObjectToInputArray, &input,
- NI_ObjectToOptionalInputArray, &labels, &indices_object))
- goto exit;
-
- if (!_NI_GetIndices(indices_object, &result_indices, &min_label,
- &max_label, &n_results))
- goto exit;
-
- center_of_mass = (double*)malloc(input->nd * n_results *
- sizeof(double));
- if (!center_of_mass) {
- PyErr_NoMemory();
- goto exit;
- }
-
- if (!NI_CenterOfMass(input, labels, min_label, max_label,
- result_indices, n_results, center_of_mass))
- goto exit;
-
- result = _NI_BuildMeasurementResultDoubleTuple(n_results, input->nd,
- center_of_mass);
-
- exit:
- Py_XDECREF(input);
- Py_XDECREF(labels);
- if (result_indices)
- free(result_indices);
- if (center_of_mass)
- free(center_of_mass);
- return result;
-}
-
-static PyObject *Py_Histogram(PyObject *obj, PyObject *args)
-{
- PyArrayObject *input = NULL, *labels = NULL, **histograms = NULL;
- PyObject *indices_object, *result = NULL;
- maybelong min_label, max_label, *result_indices = NULL, n_results;
- maybelong jj, nbins;
- long nbins_in;
- double min, max;
-
- if (!PyArg_ParseTuple(args, "O&ddlO&O", NI_ObjectToInputArray, &input,
- &min, &max, &nbins_in, NI_ObjectToOptionalInputArray,
- &labels, &indices_object))
- goto exit;
- nbins = nbins_in;
-
- if (!_NI_GetIndices(indices_object, &result_indices, &min_label,
- &max_label, &n_results))
- goto exit;
-
- /* Set all pointers to NULL, so that freeing the memory */
- /* doesn't cause problems. */
- histograms = (PyArrayObject**)calloc(input->nd * n_results,
- sizeof(PyArrayObject*));
- if (!histograms) {
- PyErr_NoMemory();
- goto exit;
- }
- for(jj = 0; jj < n_results; jj++) {
- histograms[jj] = NA_NewArray(NULL, tInt32, 1, &nbins);
- if (!histograms[jj]) {
- PyErr_NoMemory();
- goto exit;
- }
- }
-
- if (!NI_Histogram(input, labels, min_label, max_label, result_indices,
- n_results, histograms, min, max, nbins))
- goto exit;
-
- result = _NI_BuildMeasurementResultArrayObject(n_results, histograms);
-
- exit:
- Py_XDECREF(input);
- Py_XDECREF(labels);
- if (result_indices)
- free(result_indices);
- if (histograms) {
- for(jj = 0; jj < n_results; jj++) {
- Py_XDECREF(histograms[jj]);
- }
- free(histograms);
- }
- return result;
-}
-
static PyObject *Py_DistanceTransformBruteForce(PyObject *obj,
PyObject *args)
{
@@ -1293,12 +900,6 @@
METH_VARARGS, NULL},
{"watershed_ift", (PyCFunction)Py_WatershedIFT,
METH_VARARGS, NULL},
- {"statistics", (PyCFunction)Py_Statistics,
- METH_VARARGS, NULL},
- {"center_of_mass", (PyCFunction)Py_CenterOfMass,
- METH_VARARGS, NULL},
- {"histogram", (PyCFunction)Py_Histogram,
- METH_VARARGS, NULL},
{"distance_transform_bf", (PyCFunction)Py_DistanceTransformBruteForce,
METH_VARARGS, NULL},
{"distance_transform_op", (PyCFunction)Py_DistanceTransformOnePass,
Modified: trunk/scipy/ndimage/tests/test_ndimage.py
===================================================================
--- trunk/scipy/ndimage/tests/test_ndimage.py 2010-04-23 00:21:57 UTC (rev 6338)
+++ trunk/scipy/ndimage/tests/test_ndimage.py 2010-04-24 20:59:32 UTC (rev 6339)
@@ -48,6 +48,8 @@
a = numpy.asarray(a, numpy.complex128)
b = numpy.asarray(b, numpy.complex128)
t = ((a.real - b.real)**2).sum() + ((a.imag - b.imag)**2).sum()
+ if (a.dtype == numpy.object or b.dtype == numpy.object):
+ t = sum([diff(c,d)**2 for c,d in zip(a,b)])
else:
a = numpy.asarray(a)
a = a.astype(numpy.float64)
@@ -2777,15 +2779,8 @@
input = numpy.array([[1, 2], [3, 4]], type)
output = ndimage.sum(input, labels = labels,
index = [4, 8, 2])
- self.failUnless(output == [4.0, 0.0, 5.0])
+ self.failUnless(numpy.all(output == [4.0, 0.0, 5.0]))
- def test_sum13(self):
- "sum 13"
- input = numpy.array([1,2,3,4])
- labels = numpy.array([0,0,0,0])
- index = numpy.array([0],numpy.uint64)
- self.failUnlessRaises(ValueError,ndimage.sum,input,labels,index)
-
def test_mean01(self):
"mean 1"
labels = numpy.array([1, 0], bool)
@@ -2817,7 +2812,8 @@
input = numpy.array([[1, 2], [3, 4]], type)
output = ndimage.mean(input, labels = labels,
index = [4, 8, 2])
- self.failUnless(output == [4.0, 0.0, 2.5])
+ self.failUnless(numpy.all(output[[0,2]] == [4.0, 2.5]) and
+ numpy.isnan(output[1]))
def test_minimum01(self):
"minimum 1"
@@ -2850,7 +2846,7 @@
input = numpy.array([[1, 2], [3, 4]], type)
output = ndimage.minimum(input, labels = labels,
index = [2, 3, 8])
- self.failUnless(output == [2.0, 4.0, 0.0])
+ self.failUnless(numpy.all(output == [2.0, 4.0, 0.0]))
def test_maximum01(self):
"maximum 1"
@@ -2883,7 +2879,7 @@
input = numpy.array([[1, 2], [3, 4]], type)
output = ndimage.maximum(input, labels = labels,
index = [2, 3, 8])
- self.failUnless(output == [3.0, 4.0, 0.0])
+ self.failUnless(numpy.all(output == [3.0, 4.0, 0.0]))
def test_maximum05(self):
"Ticket #501"
@@ -2895,7 +2891,7 @@
for type in self.types:
input = numpy.array([], type)
output = ndimage.variance(input)
- self.failUnless(float(output) == 0.0)
+ self.failUnless(numpy.isnan(output))
def test_variance02(self):
"variance 2"
@@ -2909,13 +2905,13 @@
for type in self.types:
input = numpy.array([1, 3], type)
output = ndimage.variance(input)
- self.failUnless(output == 2.0)
+ self.failUnless(output == 1.0)
def test_variance04(self):
"variance 4"
input = numpy.array([1, 0], bool)
output = ndimage.variance(input)
- self.failUnless(output == 0.5)
+ self.failUnless(output == 0.25)
def test_variance05(self):
"variance 5"
@@ -2923,7 +2919,7 @@
for type in self.types:
input = numpy.array([1, 3, 8], type)
output = ndimage.variance(input, labels, 2)
- self.failUnless(output == 2.0)
+ self.failUnless(output == 1.0)
def test_variance06(self):
"variance 6"
@@ -2931,14 +2927,14 @@
for type in self.types:
input = numpy.array([1, 3, 8, 10, 8], type)
output = ndimage.variance(input, labels, [2, 3, 4])
- self.failUnless(output == [2.0, 2.0, 0.0])
+ self.failUnless(numpy.all(output == [1.0, 1.0, 0.0]))
def test_standard_deviation01(self):
"standard deviation 1"
for type in self.types:
input = numpy.array([], type)
output = ndimage.standard_deviation(input)
- self.failUnless(float(output) == 0.0)
+ self.failUnless(numpy.isnan(output))
def test_standard_deviation02(self):
"standard deviation 2"
@@ -2952,13 +2948,13 @@
for type in self.types:
input = numpy.array([1, 3], type)
output = ndimage.standard_deviation(input)
- self.failUnless(output == math.sqrt(2.0))
+ self.failUnless(output == math.sqrt(1.0))
def test_standard_deviation04(self):
"standard deviation 4"
input = numpy.array([1, 0], bool)
output = ndimage.standard_deviation(input)
- self.failUnless(output == math.sqrt(0.5))
+ self.failUnless(output == 0.5)
def test_standard_deviation05(self):
"standard deviation 5"
@@ -2966,7 +2962,7 @@
for type in self.types:
input = numpy.array([1, 3, 8], type)
output = ndimage.standard_deviation(input, labels, 2)
- self.failUnless(output == math.sqrt(2.0))
+ self.failUnless(output == 1.0)
def test_standard_deviation06(self):
"standard deviation 6"
@@ -2975,8 +2971,7 @@
input = numpy.array([1, 3, 8, 10, 8], type)
output = ndimage.standard_deviation(input, labels,
[2, 3, 4])
- self.failUnless(output == [math.sqrt(2.0), math.sqrt(2.0),
- 0.0])
+ self.failUnless(all(output == [1.0, 1.0, 0.0]))
def test_minimum_position01(self):
"minimum position 1"
@@ -3041,7 +3036,7 @@
[1, 5, 1, 1]], type)
output = ndimage.minimum_position(input, labels,
[2, 3])
- self.failUnless(output == [(0, 1), (1, 2)])
+ self.failUnless(output[0] == (0, 1) and output[1] == (1, 2))
def test_maximum_position01(self):
"maximum position 1"
@@ -3098,7 +3093,7 @@
[1, 5, 1, 1]], type)
output = ndimage.maximum_position(input, labels,
[1, 2])
- self.failUnless(output == [(0, 0), (1, 1)])
+ self.failUnless(output[0] == (0, 0) and output[1] == (1, 1))
def test_extrema01(self):
"extrema 1"
@@ -3148,8 +3143,10 @@
labels = labels, index = [2, 3, 8])
output5 = ndimage.maximum_position(input,
labels = labels, index = [2, 3, 8])
- self.failUnless(output1 == (output2, output3, output4,
- output5))
+ self.failUnless(numpy.all(output1[0] == output2))
+ self.failUnless(numpy.all(output1[1] == output3))
+ self.failUnless(numpy.all(output1[2] == output4))
+ self.failUnless(numpy.all(output1[3] == output5))
def test_extrema04(self):
"extrema 4"
@@ -3165,8 +3162,10 @@
[1, 2])
output5 = ndimage.maximum_position(input, labels,
[1, 2])
- self.failUnless(output1 == (output2, output3, output4,
- output5))
+ self.failUnless(numpy.all(output1[0] == output2))
+ self.failUnless(numpy.all(output1[1] == output3))
+ self.failUnless(numpy.all(output1[2] == output4))
+ self.failUnless(numpy.all(output1[3] == output5))
def test_center_of_mass01(self):
"center of mass 1"
@@ -3260,7 +3259,7 @@
def test_histogram02(self):
"histogram 2"
labels = [1, 1, 1, 1, 2, 2, 2, 2]
- true = [0, 2, 0, 1, 0]
+ true = [0, 2, 0, 1, 1]
input = numpy.array([1, 1, 3, 4, 3, 3, 3, 3])
output = ndimage.histogram(input, 0, 4, 5, labels, 1)
e = diff(true, output)
@@ -3269,7 +3268,7 @@
def test_histogram03(self):
"histogram 3"
labels = [1, 0, 1, 1, 2, 2, 2, 2]
- true1 = [0, 1, 0, 1, 0]
+ true1 = [0, 1, 0, 1, 1]
true2 = [0, 0, 0, 3, 0]
input = numpy.array([1, 1, 3, 4, 3, 5, 3, 3])
output = ndimage.histogram(input, 0, 4, 5, labels, (1,2))
More information about the Scipy-svn
mailing list