Reconstruction

.gitignore

@@ -1,3 +1,4 @@
+openpiv/docs/sphinx-docs/downloads/tutorial-part1
 *.m
 openpiv/src/process.c
 openpiv/src/lib.c

openpiv/docs/sphinx-docs/src/api_reference.rst

@@ -56,9 +56,20 @@ The ``openpiv.process`` module
     correlate_windows
     normalize_intensity
 
+The ``openpiv.lib`` module
+--------------------------------
+.. automodule:: openpiv.lib
+
+.. currentmodule:: openpiv.lib
+
+.. autosummary:: 
+    :toctree: generated/
+
+    sincinterp    
+    replace_nans
 
 
-The ``openpiv.filters`` module
+The ``openpiv.filters`` module  
 ------------------------------
 
 .. automodule:: openpiv.filters

openpiv/docs/sphinx-docs/src/generated/openpiv.lib.replace_nans.rst

@@ -0,0 +1,6 @@
+openpiv.lib.replace_nans
+========================
+
+.. currentmodule:: openpiv.lib
+
+.. autofunction:: replace_nans

openpiv/docs/sphinx-docs/src/generated/openpiv.lib.sincinterp.rst

@@ -0,0 +1,6 @@
+openpiv.lib.sincinterp
+======================
+
+.. currentmodule:: openpiv.lib
+
+.. autofunction:: sincinterp

openpiv/src/lib.pyx

@@ -7,6 +7,9 @@ cimport cython
 DTYPEf = np.float64
 ctypedef np.float64_t DTYPEf_t
 
+DTYPEi = np.int32
+ctypedef np.int32_t DTYPEi_t
+
 @cython.boundscheck(False) # turn of bounds-checking for entire function
 @cython.wraparound(False) # turn of bounds-checking for entire function
 def replace_nans( np.ndarray[DTYPEf_t, ndim=2] array, int max_iter, float tol, int kernel_size=1, str method='localmean'):
@@ -125,3 +128,77 @@ def replace_nans( np.ndarray[DTYPEf_t, ndim=2] array, int max_iter, float tol, i
                 replaced_old[l] = replaced_new[l]
 
     return filled
+
+
+def sincinterp( np.ndarray[DTYPEi_t, ndim=2] image, np.ndarray[DTYPEf_t, ndim=2] x, np.ndarray[DTYPEf_t, ndim=2] y, int kernel_size=3 ):
+    """Re-sample an image at intermediate positions between pixels.
+
+    This function uses a cardinal interpolation formula which limits 
+    the loss of information in the resampling process. It uses a limited
+    number of neighbouring pixels.
+
+
+    The new image :math:`im^+` at fractional locations :math:`x` and :math:`y` is computed as:
+
+    .. math:: 
+
+       im^+(x,y) = \sum_{i=-\mathtt{kernel\_size}}^{i=\mathtt{kernel\_size}} \sum_{j=-\mathtt{kernel\_size}}^{j=\mathtt{kernel\_size}} \mathtt{image}(i,j)  sin[\pi(i-\mathtt{x})]  sin[\pi(j-\mathtt{y})]  / \pi(i-\mathtt{x}) / \pi(j-\mathtt{y}) 
+
+
+    Parameters
+    ----------
+    image : np.ndarray, dtype np.int32
+        the image array.
+
+    x : two dimensions np.ndarray of floats
+        an array containing fractional pixel row
+        positions at which to interpolate the image
+
+    y : two dimensions np.ndarray of floats
+        an array containing fractional pixel column 
+        positions at which to interpolate the image
+
+    kernel_size : int
+        interpolation is performed over a ``(2*kernel_size+1)*(2*kernel_size+1)`` 
+        submatrix  in the neighbourhood of each interpolation point.
+
+    Returns
+    -------
+
+    im : np.ndarray, dtype np.float64
+        the interpolated value of ``image`` at the points specified 
+        by ``x`` and ``y``
+
+    """
+
+    # indices
+    cdef int i, j, I, J
+
+    # the output array
+    cdef np.ndarray[DTYPEf_t, ndim=2] r = np.zeros( [x.shape[0], x.shape[1]], dtype=DTYPEf)
+
+    # fast pi
+    cdef float pi = 3.1419
+
+    # for each point of the output array
+    for I in range(x.shape[0]):
+        for J in range(x.shape[1]):
+
+            #loop over all neighbouring grid points 
+            for i in range( int(x[I,J])-kernel_size, int(x[I,J])+kernel_size+1 ):
+                for j in range( int(y[I,J])-kernel_size, int(y[I,J])+kernel_size+1 ):
+                    # check that we are in the boundaries
+                    if i >= 0 and i <= image.shape[0] and j >= 0 and j <= image.shape[1]:
+                        if (i-x[I,J]) == 0.0 and (j-y[I,J]) == 0.0:
+                            r[I,J] = r[I,J] + image[i,j]
+                        elif (i-x[I,J]) == 0.0:
+                            r[I,J] = r[I,J] + image[i,j] * sin( pi*(j-y[I,J]) )/( pi*(j-y[I,J]) )
+                        elif (j-y[I,J]) == 0.0:
+                            r[I,J] = r[I,J] + image[i,j] * sin( pi*(i-x[I,J]) )/( pi*(i-x[I,J]) )
+                        else:
+                            r[I,J] = r[I,J] + image[i,j] * sin( pi*(i-x[I,J]) )*sin( pi*(j-y[I,J]) )/( pi*pi*(i-x[I,J])*(j-y[I,J]))
+    return r
+
+
+cdef extern from "math.h":
+    double sin(double)