add rank_downscale

gputools/convolve/__init__.py

@@ -12,3 +12,4 @@
 #from .minmax_filter import max_filter, min_filter
 from .generic_separable_filters import max_filter, min_filter, uniform_filter
 from .median_filter import median_filter
+from .rank_downscale import rank_downscale

gputools/convolve/kernels/rank_downscale.cl

@@ -0,0 +1,124 @@
+// https://github.com/rudrasohan/algorithms/blob/master/quickselect.c
+
+inline void swap(${DTYPE}* a, ${DTYPE}* b)
+{
+    ${DTYPE} t = *a;
+    *a = *b;
+    *b = t;
+}
+
+inline int partition(${DTYPE}* arr, const int l, const int r) 
+{ 
+    ${DTYPE} x = arr[r];
+    int i = l; 
+    for (int j = l; j <= r - 1; j++) { 
+        if (arr[j] <= x) { 
+            swap(&arr[i], &arr[j]); 
+            i++; 
+        } 
+    } 
+    swap(&arr[i], &arr[r]); 
+    return i; 
+} 
+
+inline ${DTYPE} kthSmallest(${DTYPE} arr[], int l, int r, int k) 
+{ 
+    // If k is smaller than number of  
+    // elements in array 
+    if (k > 0 && k <= r - l + 1) { 
+
+        // Partition the array around last  
+        // element and get position of pivot  
+        // element in sorted array 
+        int index = partition(arr, l, r); 
+
+        // If position is same as k 
+        if (index - l == k - 1) 
+            return arr[index]; 
+
+        // If position is more, recur  
+        // for left subarray 
+        if (index - l > k - 1)  
+            return kthSmallest(arr, l, index - 1, k); 
+
+        // Else recur for right subarray 
+        return kthSmallest(arr, index + 1, r,  
+                            k - index + l - 1); 
+    } 
+    else
+        return (${DTYPE})0;
+
+} 
+
+
+
+__kernel void rank_2(__global ${DTYPE} * input,
+                    __global ${DTYPE} * output,
+                    const int Nx0, const int Ny0, 
+                    const int rank){
+
+  int x = get_global_id(0);
+  int y = get_global_id(1);
+
+  int Nx = get_global_size(0);
+  int Ny = get_global_size(1);
+
+
+
+  ${DTYPE} a[${FSIZE_Y}*${FSIZE_X}];
+
+  for (int m = 0; m < ${FSIZE_Y}; ++m) {
+    for (int n = 0; n < ${FSIZE_X}; ++n) {
+
+	  int x2 = x*${FSIZE_X}+n;
+	  int y2 = y*${FSIZE_Y}+m;
+
+	  bool inside = ((x2>=0)&&(x2<Nx0)&&(y2>=0)&&(y2<Ny0));
+
+	  a[n+${FSIZE_X}*m] = inside?input[x2+y2*Nx0]:(${DTYPE})(${CVAL});
+	  }
+	}
+
+
+  output[x+y*Nx] = kthSmallest(a, 0, ${FSIZE_X}*${FSIZE_Y}-1, rank+1);
+
+}
+
+
+
+__kernel void rank_3(__global ${DTYPE} * input,
+						__global ${DTYPE} * output,
+                        const int Nx0, const int Ny0, const int Nz0,
+                        const int rank){
+
+  int x = get_global_id(0);
+  int y = get_global_id(1);
+  int z = get_global_id(2);
+
+  int Nx = get_global_size(0);
+  int Ny = get_global_size(1);
+  int Nz = get_global_size(2);
+
+
+
+  ${DTYPE} a[${FSIZE_Z}*${FSIZE_Y}*${FSIZE_X}];
+
+  for (int p = 0; p < ${FSIZE_Z}; ++p) {
+	for (int m = 0; m < ${FSIZE_Y}; ++m) {
+	  for (int n = 0; n < ${FSIZE_X}; ++n) {
+
+	  int x2 = x*${FSIZE_X}+n;
+	  int y2 = y*${FSIZE_Y}+m;
+	  int z2 = z*${FSIZE_Z}+p;
+
+	  bool inside = ((x2>=0)&&(x2<Nx0)&&(y2>=0)&&(y2<Ny0)&&(z2>=0)&&(z2<Nz0));
+
+	  a[n+${FSIZE_X}*m+${FSIZE_X}*${FSIZE_Y}*p] = inside?input[x2+y2*Nx0+z2*Nx0*Ny0]:(${DTYPE})(${CVAL});
+	  }
+	}
+  }
+
+
+  output[x+y*Nx+z*Nx*Ny] = kthSmallest(a, 0, ${FSIZE_X}*${FSIZE_Y}*${FSIZE_Z}-1, rank+1);
+
+}

gputools/convolve/rank_downscale.py

@@ -0,0 +1,114 @@
+from __future__ import print_function, unicode_literals, absolute_import, division
+import logging
+from typing import Union, Tuple
+
+logger = logging.getLogger(__name__)
+
+import numpy as np
+from mako.template import Template
+from numbers import Number
+from gputools import OCLArray, OCLProgram
+from gputools.core.ocltypes import cl_buffer_datatype_dict
+from ._abspath import abspath
+
+
+def _rank_downscale_2d(data_g, size=(3, 3), rank=None, cval = 0, res_g=None):
+    if not data_g.dtype.type in cl_buffer_datatype_dict:
+        raise ValueError("dtype %s not supported" % data_g.dtype.type)
+
+    DTYPE = cl_buffer_datatype_dict[data_g.dtype.type]
+
+    size = tuple(map(int, size)) 
+
+    if rank is None:
+        rank = np.prod(size) // 2
+
+    with open(abspath("kernels/rank_downscale.cl"), "r") as f:
+        tpl = Template(f.read())
+
+    rendered = tpl.render(DTYPE = DTYPE,FSIZE_Z=0, FSIZE_X=size[1], FSIZE_Y=size[0],CVAL = cval)
+
+    prog = OCLProgram(src_str=rendered)
+
+    if res_g is None:
+        res_g = OCLArray.empty(tuple(s0//s for s, s0 in zip(size,data_g.shape)), data_g.dtype)
+
+    prog.run_kernel("rank_2", res_g.shape[::-1], None, data_g.data, res_g.data, 
+                    np.int32(data_g.shape[1]), np.int32(data_g.shape[0]),
+                    np.int32(rank))
+    return res_g
+
+def _rank_downscale_3d(data_g, size=(3, 3, 3), rank=None, cval = 0, res_g=None):
+    if not data_g.dtype.type in cl_buffer_datatype_dict:
+        raise ValueError("dtype %s not supported" % data_g.dtype.type)
+
+    DTYPE = cl_buffer_datatype_dict[data_g.dtype.type]
+
+    size = tuple(map(int, size)) 
+
+    if rank is None:
+        rank = np.prod(size) // 2
+
+    with open(abspath("kernels/rank_downscale.cl"), "r") as f:
+        tpl = Template(f.read())
+
+    rendered = tpl.render(DTYPE = DTYPE,FSIZE_X=size[2], FSIZE_Y=size[1], FSIZE_Z=size[0],CVAL = cval)
+
+    prog = OCLProgram(src_str=rendered)
+
+    if res_g is None:
+        res_g = OCLArray.empty(tuple(s0//s for s, s0 in zip(size,data_g.shape)), data_g.dtype)
+
+    prog.run_kernel("rank_3", res_g.shape[::-1], None, data_g.data, res_g.data, 
+                    np.int32(data_g.shape[2]), np.int32(data_g.shape[1]), np.int32(data_g.shape[0]),
+                    np.int32(rank))
+    return res_g
+
+
+def rank_downscale(data:np.ndarray, size:Union[int, Tuple[int]]=3, rank=None):
+    """
+        downscales an image by the given factor and returns the rank-th element in each block
+
+    Parameters
+    ----------
+    data: numpy.ndarray
+        input data (2d or 3d)
+    size: int or tuple
+        downsampling factors
+    rank: int
+        rank of element to retain 
+            rank = 0     -> minimum 
+            rank = -1    -> maximum 
+            rank = None  -> median 
+        
+    Returns
+    -------
+        downscaled image 
+    """
+
+    if not (isinstance(data, np.ndarray) and data.ndim in (2,3)):
+        raise ValueError("input data has to be a 2d or 3d numpy array!")
+
+    if isinstance(size, Number):
+        size = (int(size),)*data.ndim
+
+    if len(size) != data.ndim:
+        raise ValueError("size has to be a tuple of 3 ints")
+
+    if rank is None:
+        rank = np.prod(size) // 2
+    else:
+        rank = rank % np.prod(size)
+
+    data_g = OCLArray.from_array(data)
+
+    if data.ndim==2:
+        res_g = _rank_downscale_2d(data_g, size=size, rank=rank)
+    elif data.ndim==3:
+        res_g = _rank_downscale_3d(data_g, size=size, rank=rank)
+    else: 
+        raise ValueError("data has to be 2d or 3d")
+
+    return res_g.get()
+
+

gputools/version.py

@@ -6,4 +6,4 @@
 
 """
 
-__version__ = "0.2.15"
+__version__ = "0.2.16"