Skip to content
This repository

ENH: add sum_angle() and sum_polar() functions to twodim_base.py #230

Open
wants to merge 1 commit into from

5 participants

Robert Jordens Ralf Gommers Charles Harris Travis E. Oliphant Stefan van der Walt
Robert Jordens

sum_angle() computes the sum of a 2-d array along an angled axis
sum_polar() along radial lines or along azimuthal circles

Robert Jordens add sum_angle() and sum_polar() functions.
sum_angle() computes the sum of a 2-d array along an angled axis
sum_polar() along radial lines or along azimuthal circles
ba7766e
Ralf Gommers
Owner

Hi @jordens, I can see how this would be useful, but adding new functions should really be discussed on the numpy-discussion mailing list first. Could you please bring it up there?

Charles Harris charris commented on the diff May 20, 2012
numpy/lib/twodim_base.py
((4 lines not shown))
  885
+
  886
+
  887
+def sum_angle(m, angle, aspect=1., binsize=None):
  888
+    """
  889
+    Compute the sum of a 2-d array along an rotated axis.
  890
+
  891
+    Parameters
  892
+    ----------
  893
+    m : array_like, shape(N, M)
  894
+        2-d input array to be summed
  895
+    angle : float
  896
+        The angle of the summation direction defined such that:
  897
+            ``sum_angle(m, angle=0) == sum(m, axis=0)``
  898
+            ``sum_angle(m, angle=pi/2) == sum(m, axis=1)``
  899
+    aspect : float, optional
  900
+        The input bin aspect ratio (second dimension/first dimension)
2
Charles Harris Owner
charris added a note May 20, 2012

Dimension or component?

Robert Jordens
jordens added a note May 20, 2012

Dimension. Or I don`t understand what you mean by component.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Charles Harris
Owner
charris commented May 20, 2012

Correct me if I am wrong, but this looks like it is intended to work on images. Is that the case?

Robert Jordens
jordens commented May 20, 2012

Images are certainly a major application. But this is also interesting in the analysis of matrices that represent graphs.

Charles Harris
Owner
charris commented May 20, 2012

I was asking because I thought this might be more appropriate for one of the image processing packages. In particular, it seems related to the Hough transform. I assume the application to graphs uses the adjacency matrix. In any case, it might be worth putting up a post on the scipy mailing list as well as here since that is where graph and image algorithms are located.

Travis E. Oliphant
Owner

The sum_angle function seems equivalent to the scipy.misc.radon function which is deprecated (but a very useful function). Because your implementation is so straightforward and well documented it could be useful in NumPy as well, but I agree with Chuck that you should ask on the SciPy list as well.

Robert Jordens
jordens commented May 20, 2012

Yes. angle_sum is a Hough or Radon transform. The implementation in scipy has a couple of problems: It is based on PIL and does not seem to work with floats. It does not conserve the sum over all entries (independent of padding...). Finally, the interpolations used in PIL do not really make sense in this application. It is also 40% slower than this one.

Scikits-image has taken the radon transform from scipy: https://github.com/scikits-image/scikits-image/blob/master/skimage/transform/radon_transform.py I believe that one does not suffer the problems the PIL-based one has. But it looks even slower.

I'll raise this on the scipy/scikits list.

Stefan van der Walt
Collaborator
stefanv commented May 21, 2012

@jordens Could we use your sum_angle in scikits-image? I'd love to improve the execution speed of the radon and hough transforms (actually, I think the hough tf already uses the bincount trick, so it may not benefit). Perhaps you would be interested in working on a PR together?

Robert Jordens
jordens commented May 21, 2012

@stefanv Sure. Feel free. I will help if I can find the time.

OTOH I would still love to see these added to numpy and agree with Travis' reasoning.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Showing 1 unique commit by 1 author.

Mar 07, 2012
Robert Jordens add sum_angle() and sum_polar() functions.
sum_angle() computes the sum of a 2-d array along an angled axis
sum_polar() along radial lines or along azimuthal circles
ba7766e
This page is out of date. Refresh to see the latest.

Showing 1 changed file with 231 additions and 2 deletions. Show diff stats Hide diff stats

  1. 233  numpy/lib/twodim_base.py
233  numpy/lib/twodim_base.py
@@ -5,11 +5,12 @@
5 5
 __all__ = ['diag','diagflat','eye','fliplr','flipud','rot90','tri','triu',
6 6
            'tril','vander','histogram2d','mask_indices',
7 7
            'tril_indices','tril_indices_from','triu_indices','triu_indices_from',
8  
-           ]
  8
+           'sum_angle', 'sum_polar']
9 9
 
10 10
 from numpy.core.numeric import asanyarray, equal, subtract, arange, \
11 11
      zeros, greater_equal, multiply, ones, asarray, alltrue, where, \
12  
-     empty, diagonal
  12
+     empty, diagonal, sin, cos, arctan2, pi, floor
  13
+from numpy.lib.function_base import bincount
13 14
 
14 15
 def fliplr(m):
15 16
     """
@@ -881,3 +882,231 @@ def triu_indices_from(arr, k=0):
881 882
     if not (arr.ndim == 2 and arr.shape[0] == arr.shape[1]):
882 883
         raise ValueError("input array must be 2-d and square")
883 884
     return triu_indices(arr.shape[0],k)
  885
+
  886
+
  887
+def sum_angle(m, angle, aspect=1., binsize=None):
  888
+    """
  889
+    Compute the sum of a 2-d array along an rotated axis.
  890
+
  891
+    Parameters
  892
+    ----------
  893
+    m : array_like, shape(N, M)
  894
+        2-d input array to be summed
  895
+    angle : float
  896
+        The angle of the summation direction defined such that:
  897
+            ``sum_angle(m, angle=0) == sum(m, axis=0)``
  898
+            ``sum_angle(m, angle=pi/2) == sum(m, axis=1)``
  899
+    aspect : float, optional
  900
+        The input bin aspect ratio (second dimension/first dimension)
  901
+    binsize : float, optional
  902
+        The output bin size in units of the first input dimension step
  903
+        size. If no binsize is given, it defaults to the "natural bin
  904
+        size" which is the larger projection of the two input step sizes
  905
+        onto the output dimension (the axis perpendicular to the
  906
+        summation axis).
  907
+
  908
+    Returns
  909
+    -------
  910
+    out : ndarray, shape(K)
  911
+        The sum of `m` along the axis at `angle`
  912
+
  913
+    See Also
  914
+    --------
  915
+    sum_polar : similar method summing azimuthally or radially
  916
+
  917
+    Notes
  918
+    -----
  919
+    The summation angle is relative to the first dimension.
  920
+
  921
+    For 0<=angle<=pi/2 the value at [0,0] ends up in the first bin and
  922
+    the value at [-1,-1] ends up in the last bin. Up to rounding, the
  923
+    center value will always end up in the center bin.
  924
+
  925
+    For angle=3/4*pi the summation is along the diagonal.
  926
+    For angle=1/4*pi the summation is along the antidiagonal.
  927
+
  928
+    The origin of the rotation is the [0,0] index. This determines the
  929
+    bin rounding.
  930
+
  931
+    There is no interpolation and artefacts are likely if this function
  932
+    is interpreted as an image processing function.
  933
+
  934
+    The full array sum is always strictly conserved:
  935
+        ``sum_angle(m, t).sum() == m.sum()``
  936
+
  937
+    .. versionadded:: 1.4.0
  938
+
  939
+    Examples
  940
+    --------
  941
+    >>> m = np.arange(9.).reshape((3, 3))
  942
+    >>> np.all(sum_angle(m, 0) == m.sum(axis=0))
  943
+    True
  944
+    >>> np.all(sum_angle(m, np.pi/2) == m.sum(axis=1))
  945
+    True
  946
+    >>> np.all(sum_angle(m, np.pi) == m.sum(axis=0)[::-1])
  947
+    True
  948
+    >>> np.all(sum_angle(m, 3*np.pi/2) == m.sum(axis=1)[::-1])
  949
+    True
  950
+    >>> np.all(sum_angle(m, 2*np.pi) == m.sum(axis=0))
  951
+    True
  952
+    >>> np.all(sum_angle(m, -np.pi/2) ==
  953
+    ...        sum_angle(m, 3*np.pi/2))
  954
+    True
  955
+    >>> d1 = np.array([0, 4, 12, 12, 8]) # antidiagonal
  956
+    >>> d2 = np.array([2, 6, 12, 10, 6]) # diagonal
  957
+    >>> np.all(sum_angle(m, np.pi/4) == d1)
  958
+    True
  959
+    >>> np.all(sum_angle(m, 3*np.pi/4) == d2)
  960
+    True
  961
+    >>> np.all(sum_angle(m, 5*np.pi/4) == d1[::-1])
  962
+    True
  963
+    >>> np.all(sum_angle(m, 7*np.pi/4) == d2[::-1])
  964
+    True
  965
+    >>> np.all(sum_angle(m, 0, aspect=2, binsize=1) ==
  966
+    ...        np.array([9, 0, 12, 0, 15]))
  967
+    True
  968
+    >>> np.all(sum_angle(m, 0, aspect=.5, binsize=1) ==
  969
+    ...        np.array([9, 12+15]))
  970
+    True
  971
+    >>> np.all(sum_angle(m, 0, aspect=.5) == m.sum(axis=0))
  972
+    True
  973
+    >>> np.all(sum_angle(m, np.pi/2, aspect=2, binsize=1) ==
  974
+    ...        m.sum(axis=1))
  975
+    True
  976
+    >>> m2 = np.arange(1e6).reshape((100, 10000))
  977
+    >>> np.all(sum_angle(m2, 0) == m2.sum(axis=0))
  978
+    True
  979
+    >>> np.all(sum_angle(m2, np.pi/2) == m2.sum(axis=1))
  980
+    True
  981
+    >>> sum_angle(m2, np.pi/4).shape
  982
+    (10099,)
  983
+    >>> sum_angle(m2, np.pi/4).sum() == m2.sum()
  984
+    True
  985
+    """
  986
+    m = asanyarray(m)
  987
+    if m.ndim != 2:
  988
+        raise ValueError("Input must be 2-d.")
  989
+    if binsize is None:
  990
+        binsize = max(abs(cos(angle)*aspect), abs(sin(angle)))
  991
+    # first axis needs to be inverted due to the angle convention
  992
+    m = m[::-1]
  993
+    i, j = arange(m.shape[0])[:, None], np.arange(m.shape[1])[None, :]
  994
+    k = (cos(angle)*aspect/binsize)*j - (sin(angle)/binsize)*i
  995
+    cx, cy = (0, 0, -1, -1), (0, -1, 0, -1)
  996
+    km = k[cx, cy].min()
  997
+    # kp = k[cx, cy].max()
  998
+    # minlength=kp-km
  999
+    k = floor(k-(km-.5)).astype(int)
  1000
+    return bincount(k.ravel(), m.ravel())
  1001
+
  1002
+
  1003
+def sum_polar(m, center, direction, aspect=1., binsize=None):
  1004
+    """
  1005
+    Compute the sum of a 2-d array radially or azimuthally.
  1006
+
  1007
+    Parameters
  1008
+    ----------
  1009
+    m : array_like, shape(N, M)
  1010
+        2-d input array to be summed
  1011
+    center : tuple(float, float)
  1012
+        The center of the summation measured from the [0, 0] index
  1013
+        in units of the two input step sizes
  1014
+    direction : "radial" or "azimuthal"
  1015
+        Summation direction
  1016
+    aspect : float, optional
  1017
+        The input bin aspect ratio (second dimension/first dimension)
  1018
+    binsize : int, optional
  1019
+        The output bin size. If None is given, and direction="radial"
  1020
+        then binsize=2*pi/100, else binsize=min(1, aspect).
  1021
+
  1022
+    Returns
  1023
+    -------
  1024
+    out : ndarray, shape(K)
  1025
+        The radial or azimuthal sum of `m`
  1026
+
  1027
+    See Also
  1028
+    --------
  1029
+    sum_angle : similar method summing along angled parallel lines
  1030
+
  1031
+    Notes
  1032
+    -----
  1033
+    The index of `out` is the binned radius or the binned angle, both
  1034
+    according to `binsize`.
  1035
+
  1036
+    Angles are measured from the positive second index axis towards the
  1037
+    negative first index axis. They correspond to mathematically
  1038
+    positive angles in the index coordinates of m[::-1] -- or the [0, 0]
  1039
+    index in the lower left.
  1040
+
  1041
+    If direction="azimuthal" then the length of the output is determined
  1042
+    by the maximal distance to the center. The radius-bins are [0,
  1043
+    binsize), [binsize, 2*binsize), ... up to [r, r+binsize) for some
  1044
+    value r with max_radius-binsize <= r < max_radius.
  1045
+
  1046
+    If direction="radial" the length is always 2*pi/binsize. This is not
  1047
+    the same as arctan2(i, j) which would distinguish +pi and -pi!  The
  1048
+    azimuthal bins are therefore [-pi, -pi+binsize), [-pi+binsize,
  1049
+    2*binsize), ... up to [p-binsize, p) for some p with pi-binsize <= p
  1050
+    < pi. The values at +pi end up in the first bin.  See ``arctan2``
  1051
+    for the behaviour in other special cases.
  1052
+
  1053
+    There is no interpolation and artefacts are likely if this function
  1054
+    is interpreted as an image processing function.
  1055
+
  1056
+    The full array sum is always strictly conserved:
  1057
+        ``sum_polar(m, ...).sum() == m.sum()``
  1058
+
  1059
+    .. versionadded:: 1.4.0
  1060
+
  1061
+    Examples
  1062
+    --------
  1063
+    >>> m = np.arange(1., 10.).reshape((3, 3))
  1064
+    >>> sum_polar(m, (0, 0), "radial").sum() == m.sum()
  1065
+    True
  1066
+    >>> sum_polar(m, (0, 0), "azimuthal").sum() == m.sum()
  1067
+    True
  1068
+    >>> sum_polar(m, (1, 1), "radial").sum() == m.sum()
  1069
+    True
  1070
+    >>> sum_polar(m, (1, 1), "azimuthal").sum() == m.sum()
  1071
+    True
  1072
+    >>> sum_polar(m, (1, 1), "radial", binsize=np.pi/4)
  1073
+    array([  4.,   1.,   2.,   3.,  11.,   9.,   8.,   7.])
  1074
+    >>> sum_polar(m, (1, 1), "azimuthal", binsize=1.)
  1075
+    array([  5.,  40.])
  1076
+    >>> sum_polar(m, (1, 1), "azimuthal", binsize=2**.5/2)
  1077
+    array([  5.,  20.,  20.])
  1078
+    >>> sum_polar(m, (.5, .5), "azimuthal", binsize=1)
  1079
+    array([ 12.,  24.,   9.])
  1080
+    >>> sum_polar(m, (0, 0), "radial", binsize=np.pi/8)
  1081
+    array([  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   6.,   6.,  22.,
  1082
+             0.,  11.,   0.,   0.,   0.])
  1083
+    >>> sum_polar(m, (0, 0), "radial", binsize=np.pi/2)
  1084
+    array([  0.,   0.,  34.,  11.])
  1085
+    >>> m2 = np.arange(123*345).reshape((123, 345))
  1086
+    >>> sum_polar(m2, (67, 89), "radial", binsize=2*np.pi/1011).shape[0]
  1087
+    1011
  1088
+    """
  1089
+    m = asanyarray(m)
  1090
+    if m.ndim != 2:
  1091
+        raise ValueError("Input must be 2-d.")
  1092
+    i, j = arange(m.shape[0])[:, None], np.arange(m.shape[1])[None, :]
  1093
+    i, j = i-center[0], j-center[1]
  1094
+    if direction == "azimuthal":
  1095
+        k = (j**2*aspect**2+i**2)**.5
  1096
+        if binsize is None:
  1097
+            binsize = min(1., aspect)
  1098
+        minlength = None
  1099
+    elif direction == "radial":
  1100
+        k = arctan2(i, j*aspect)+pi
  1101
+        if binsize is None:
  1102
+            binsize = 2*pi/100
  1103
+        minlength = int(2*pi/binsize)+1
  1104
+    else:
  1105
+        raise ValueError("direction needs to be 'radial' or 'azimuthal'")
  1106
+    k = (k/binsize).astype(int)
  1107
+    r = bincount(k.ravel(), m.ravel(), minlength)
  1108
+    if direction == "radial":
  1109
+        assert r.shape[0] == minlength, (r.shape, minlength)
  1110
+        r[0] += r[-1]
  1111
+        r = r[:-1]
  1112
+    return r
Commit_comment_tip

Tip: You can add notes to lines in a file. Hover to the left of a line to make a note

Something went wrong with that request. Please try again.