Python scipy.ndimage.spline_filter() Examples

def save_spline_arrays(self, pickleName, what):
        if _DEBUG:
            print('. Pickling splines to\n{0}'.format(pickleName))
        splines = []
        for ia, a in enumerate(what):
            a = np.concatenate((a[:, self.extraRows:0:-1, :], a), axis=1)
            a = np.concatenate((a[:, :, self.extraRows:0:-1], a), axis=2)
            if self.order == 3:
                spline = ndimage.spline_filter(a)
            else:
                spline = a
            splines.append(spline)
        Imax = np.max(what[0])
        with open(pickleName, 'wb') as f:
            pickle.dump((Imax, splines), f, protocol=2)
        return splines, Imax 

def test_affine_transform09(self):
        data = numpy.array([[4, 1, 3, 2],
                               [7, 6, 8, 5],
                               [3, 5, 3, 6]])
        for order in range(0, 6):
            if (order > 1):
                filtered = ndimage.spline_filter(data,
                                                           order=order)
            else:
                filtered = data
            out = ndimage.affine_transform(filtered,[[1, 0],
                                                               [0, 1]],
                                  [-1, -1], order=order, prefilter=False)
            assert_array_almost_equal(out, [[0, 0, 0, 0],
                                       [0, 4, 1, 3],
                                       [0, 7, 6, 8]]) 

def test_affine_transform26(self):
        # test homogeneous coordinates
        data = numpy.array([[4, 1, 3, 2],
                            [7, 6, 8, 5],
                            [3, 5, 3, 6]])
        for order in range(0, 6):
            if (order > 1):
                filtered = ndimage.spline_filter(data, order=order)
            else:
                filtered = data
            tform_original = numpy.eye(2)
            offset_original = -numpy.ones((2, 1))
            tform_h1 = numpy.hstack((tform_original, offset_original))
            tform_h2 = numpy.vstack((tform_h1, [[0, 0, 1]]))
            out1 = ndimage.affine_transform(filtered, tform_original,
                                            offset_original.ravel(),
                                            order=order, prefilter=False)
            out2 = ndimage.affine_transform(filtered, tform_h1, order=order,
                                            prefilter=False)
            out3 = ndimage.affine_transform(filtered, tform_h2, order=order,
                                            prefilter=False)
            for out in [out1, out2, out3]:
                assert_array_almost_equal(out, [[0, 0, 0, 0],
                                                [0, 4, 1, 3],
                                                [0, 7, 6, 8]]) 

def test_spline04(self):
        for type_ in self.types:
            data = numpy.ones([4], type_)
            for order in range(2, 6):
                out = ndimage.spline_filter(data, order)
                assert_array_almost_equal(out, [1, 1, 1, 1]) 

def plot_NOM_3D(fname):
    from mpl_toolkits.mplot3d import Axes3D
    from matplotlib import cm
    from matplotlib.ticker import LinearLocator, FormatStrFormatter

    xL, yL, zL = np.loadtxt(fname+'.dat', unpack=True)
    nX = (yL == yL[0]).sum()
    nY = (xL == xL[0]).sum()
    x = xL.reshape((nY, nX))
    y = yL.reshape((nY, nX))
    z = zL.reshape((nY, nX))
    x1D = xL[:nX]
    y1D = yL[::nX]
#    z += z[::-1, :]
    zmax = abs(z).max()

    fig = plt.figure()
    ax = fig.gca(projection='3d')
    surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap=cm.coolwarm,
                           linewidth=0, antialiased=False, alpha=0.5)
    ax.set_zlim(-zmax, zmax)
    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

    fig.colorbar(surf, shrink=0.5, aspect=5)

    splineZ = ndimage.spline_filter(z.T)
    nrays = 1e3
    xnew = np.random.uniform(x1D[0], x1D[-1], nrays)
    ynew = np.random.uniform(y1D[0], y1D[-1], nrays)
    coords = np.array([(xnew-x1D[0]) / (x1D[-1]-x1D[0]) * (nX-1),
                       (ynew-y1D[0]) / (y1D[-1]-y1D[0]) * (nY-1)])
    znew = ndimage.map_coordinates(splineZ, coords, prefilter=True)
    ax.scatter(xnew, ynew, znew, c=znew, marker='o', color='gray', s=50,
               cmap=cm.coolwarm)

    fig.savefig(fname+'_3d.png')
    plt.show() 

def __init__(self, *args, **kwargs):
        roe.ToroidMirror.__init__(self, *args, **kwargs)
### here you specify the bump and its mesh ###
        if get_distorted_surface is None:
            self.limPhysX = [-5, 5]
            self.limPhysY = [-125, 125]
            self.get_surface_limits()
            return
        self.warpX, self.warpY, self.warpZ, self.distortedSurfaceName =\
            get_distorted_surface()
#        print('xyz sizes:')
#        print(self.warpX.min(), self.warpX.max())
#        print(self.warpY.min(), self.warpY.max())
#        print(self.warpZ.min(), self.warpZ.max())
        self.warpNX, self.warpNY = len(self.warpX), len(self.warpY)
        self.limPhysX = np.min(self.warpX), np.max(self.warpX)
        self.limPhysY = np.min(self.warpY), np.max(self.warpY)
        self.get_surface_limits()
        self.warpA, self.warpB = np.gradient(self.warpZ)
        dx = self.warpX[1] - self.warpX[0]
        dy = self.warpY[1] - self.warpY[0]
        self.warpA = np.arctan(self.warpA/dx)
        self.warpB = np.arctan(self.warpB/dy)
#        print(self.warpZ.shape)
#end# here you specify the bump and its mesh ###
        self.warpSplineZ = ndimage.spline_filter(self.warpZ)
        self.warpSplineA = ndimage.spline_filter(self.warpA)
        self.warpSplineB = ndimage.spline_filter(self.warpB) 

def pitch_time_deformation_spec(spec, stretch_rate=None, shift_rate=None):
    """ Based on https://github.com/f0k/ismir2018/blob/master/experiments/augment.py """
    shape = spec.shape
    random = (np.random.rand(2) - .5) * 2

    warnings.filterwarnings("ignore", "The behaviour of affine_transform with a "
                            "one-dimensional array supplied for the matrix parameter "
                            "has changed in scipy 0.18.0.", module="scipy.ndimage")

    if stretch_rate is None:
        stretch_rate = 1 + random[0] * config.MAX_STRETCHING

    if shift_rate is None:
        shift_rate = 1 + random[1] * config.MAX_SHIFTING_PITCH

    new_length = int(shape[1] * stretch_rate)

    # We can do shifting/stretching and cropping in a single affine
    # transform (including setting the upper bands to zero if we shift
    # down the signal so far that we exceed the input's nyquist rate)

    spec = spline_filter(spec.T, 2).astype(spec.dtype)
    spec = affine_transform(spec, (1 / stretch_rate, 1 / shift_rate),
                            output_shape=(new_length, shape[0]),
                            offset=(0, 0), mode='constant', order=2,
                            prefilter=False)

    # clip possible negative values introduced by the interpolation
    # delete the last frame if it is empty

    if (spec[-1] == 0).all():
        spec = spec[:-1]
    return np.maximum(spec.T, 0), stretch_rate 

def test_shift09(self):
        data = numpy.array([[4, 1, 3, 2],
                            [7, 6, 8, 5],
                            [3, 5, 3, 6]])
        for order in range(0, 6):
            if (order > 1):
                filtered = ndimage.spline_filter(data, order=order)
            else:
                filtered = data
            out = ndimage.shift(filtered, [1, 1], order=order, prefilter=False)
            assert_array_almost_equal(out, [[0, 0, 0, 0],
                                            [0, 4, 1, 3],
                                            [0, 7, 6, 8]]) 

def test_affine_transform09(self):
        data = numpy.array([[4, 1, 3, 2],
                            [7, 6, 8, 5],
                            [3, 5, 3, 6]])
        for order in range(0, 6):
            if (order > 1):
                filtered = ndimage.spline_filter(data, order=order)
            else:
                filtered = data
            out = ndimage.affine_transform(filtered, [[1, 0], [0, 1]],
                                           [-1, -1], order=order,
                                           prefilter=False)
            assert_array_almost_equal(out, [[0, 0, 0, 0],
                                            [0, 4, 1, 3],
                                            [0, 7, 6, 8]]) 

def test_spline05(self):
        for type_ in self.types:
            data = numpy.ones([4, 4], type_)
            for order in range(2, 6):
                out = ndimage.spline_filter(data, order=order)
                assert_array_almost_equal(out, [[1, 1, 1, 1],
                                                [1, 1, 1, 1],
                                                [1, 1, 1, 1],
                                                [1, 1, 1, 1]]) 

def test_spline01(self):
        for type in self.types:
            data = numpy.ones([], type)
            for order in range(2, 6):
                out = ndimage.spline_filter(data, order=order)
                assert_array_almost_equal(out, 1) 

def test_spline03(self):
        for type_ in self.types:
            data = numpy.ones([], type_)
            for order in range(2, 6):
                out = ndimage.spline_filter(data, order,
                                            output=type_)
                assert_array_almost_equal(out, 1) 

def test_spline02(self):
        for type_ in self.types:
            data = numpy.array([1], type_)
            for order in range(2, 6):
                out = ndimage.spline_filter(data, order=order)
                assert_array_almost_equal(out, [1]) 

def test_spline01(self):
        for type_ in self.types:
            data = numpy.ones([], type_)
            for order in range(2, 6):
                out = ndimage.spline_filter(data, order=order)
                assert_array_almost_equal(out, 1) 

def test_shift09(self):
        data = numpy.array([[4, 1, 3, 2],
                               [7, 6, 8, 5],
                               [3, 5, 3, 6]])
        for order in range(0, 6):
            if (order > 1):
                filtered = ndimage.spline_filter(data,
                                                           order=order)
            else:
                filtered = data
            out = ndimage.shift(filtered, [1, 1], order=order,
                                          prefilter=False)
            assert_array_almost_equal(out, [[0, 0, 0, 0],
                                       [0, 4, 1, 3],
                                       [0, 7, 6, 8]]) 

def test_spline05(self):
        for type in self.types:
            data = numpy.ones([4, 4], type)
            for order in range(2, 6):
                out = ndimage.spline_filter(data, order=order)
                assert_array_almost_equal(out, [[1, 1, 1, 1],
                                           [1, 1, 1, 1],
                                           [1, 1, 1, 1],
                                           [1, 1, 1, 1]]) 

def test_spline04(self):
        for type in self.types:
            data = numpy.ones([4], type)
            for order in range(2, 6):
                out = ndimage.spline_filter(data, order)
                assert_array_almost_equal(out, [1, 1, 1, 1]) 

def test_spline03(self):
        for type in self.types:
            data = numpy.ones([], type)
            for order in range(2, 6):
                out = ndimage.spline_filter(data, order,
                                            output=type)
                assert_array_almost_equal(out, 1) 

def test_spline02(self):
        for type in self.types:
            data = numpy.array([1])
            for order in range(2, 6):
                out = ndimage.spline_filter(data, order=order)
                assert_array_almost_equal(out, [1]) 

def correlate_frames(node_pos, mesh, img, ref, settings):
    """
    Parameters
    ----------
    node_pos : ndarray
       The position of the nodes
    mesh : Mesh
       The mesh object
    img : ndarray
       2d array containing the image frame
    ref : Reference
       The reference object
    settings : DICInput
       The settings which will be used during the analysis
   Returns
   -------
   updated node positions, current pixel values
    """

    logger = logging.getLogger(__name__)

    node_pos = np.copy(node_pos).astype(settings.precision)

    # Declare empty arrays
    pixel_pos = np.zeros((2, ref.n_pixels), dtype=settings.precision)
    dnod_x = np.zeros(mesh.n_nodes * 2)

    image_filtered = nd.spline_filter(img, order=settings.interpolation_order).transpose()

    for it in range(settings.maxit):

        # Find nodal positions within ROI
        np.dot(node_pos, ref.Nref_stack, out=pixel_pos)

        # Find pixel values for current coordinates
        Ic = nd.map_coordinates(image_filtered, pixel_pos, order=settings.interpolation_order, prefilter=False)

        # Calculate position increment as (B^T B)^-1 * (B^T*dIk) "Least squares solution"
        dnod = np.dot(ref.K, ref.I0_stack - Ic)

        # Add increment to nodal positions
        node_pos[0, :] += dnod[:mesh.n_nodes]
        node_pos[1, :] += dnod[mesh.n_nodes:]

        dnod_x += dnod

        # Check for convergence
        if np.max(np.abs(dnod)) < settings.tol:
            logger.info('Frame converged in %s iterations', it)
            return np.array((dnod_x[:mesh.n_nodes], dnod_x[mesh.n_nodes:])), Ic, True

        # Reset array values
    logger.info("Frame did not converge. Largest increment was %f pixels" % np.max(dnod))
    return np.array((dnod_x[:mesh.n_nodes], dnod_x[mesh.n_nodes:])), Ic, False