Newest build, version 0.0.1.8

IMED.egg-info/PKG-INFO

@@ -0,0 +1,85 @@
+Metadata-Version: 2.1
+Name: IMED
+Version: 0.0.1.8
+Summary: Image Euclidean Distance and Standardizing Transforms
+Home-page: https://github.com/jfelding/IMED
+Author: Jacob Felding
+Author-email: jfelding@gmail.com
+License: UNKNOWN
+Project-URL: Bug Tracker, https://github.com/jfelding/IMED/issues
+Platform: UNKNOWN
+Classifier: Programming Language :: Python :: 3
+Classifier: License :: OSI Approved :: GNU General Public License v3 (GPLv3)
+Classifier: Operating System :: OS Independent
+Requires-Python: >=3.6
+Description-Content-Type: text/markdown
+License-File: LICENSE
+
+# IMage Euclidean Distance (IMED)
+
+Matrix and Frequency implementation based on https://link.springer.com/content/pdf/10.1186/s40535-015-0014-6.pdf
+
+The Image Euclidean Distance is found by transforming images or other continous data using a convolution operation, then taking the standard pixel-wise Euclidean distance between the image. This transform is referred to as the 'Standardizing Transform', or ST. 
+
+This package contains five implementations of the standardizing transform. Two of them ('full', 'sep') are matrix methods that apply linear convolution, while the remainder are frequency/DFT methods ('fft', 'dct', 'dct_by_fft'). the 'fft' performs the transform using circular convolution, while 'dct' and 'dct_by_fft' give identical results and apply symmetric convolution.
+
+The most natural boundaries are often obtained using 'dct' and 'dct_by_fft' methods. The linear convolution methods tend to underestimate the boundaries, as they correspond to zero-padding the image berfore applying a certain Gaussian filter. 'fft' tends to give periodic boundary effects.
+
+The frequency methods apply an n-dimensional transform and so can be used for continous signals with any number of dimensions. The matrix methods are 2D only.
+
+The 'full' method is the original method, and is only here for completeness. It consumes a lot of memory, and is very slow. Its use is not recommended.
+
+## Getting Started with the Standardizing Transform
+**Installation**:
+Install the latest release from pypi:
+
+    pip install IMED
+
+To get started, IMED.ST_all contains a wrapper function standardizingTrans(imgs,sigma,method,eps=0,inverse=False)
+Here is the doc:
+
+    Takes n-dimensional data and returns the n-dimensional Standardized Transform.
+    Methods 'full' and 'sep' are 2D methods only.
+
+    Parameters:
+    * imgs is a signal 
+    * sigma (float)/array-like determines the zero-mean Gaussian that defines the IMED matrix G - not G^(1/2).
+      If sigma is array-like it should contain the same number of values as the number of dimensions of imgs.
+
+    * eps (float) is an optional small parameter to offset the Gaussian so that it is always numerically non-zero. 
+    This can allow deconvolution without significant noise amplification.
+    * method (string) is the method used to perform the standardizing transform. Choose between:
+     1. **'full':** Full Dense Matrix $z_{ST}= G^{1/2}z$ using eigenvalue decomposition
+     2. **'sep'** Separated Dense Matrices $z_{ST}= G_x^{1/2}z G_y^{1/2}$ using eigenvalue decomposition 
+     3. **'fft'**: Performs circular convolution using discrete fourier transforms of image and Gaussian 
+     without enforcing symmetric boundary conditions
+     4. **'dct_by_fft'**: Performs circular convolution using discrete fourier transforms of mirrored image and 
+     Gaussian to ensure symmetric boundary conditions and reduce edge effects from 'Gibbs-like phenomenon'
+     5. **'dct'**: Performs symmetric convolution using discrete cosine transform of image, which is identical to
+     6. the 'dct_by_fft  method, but should be more efficient
+    """
+
+## Parallelization of frequency methods and backend change
+The frequency methods of the standardizing transform ('fft', 'dct' and 'dct_by_fft') are implemented using scipy.fft.
+By chopping up the transforms into smaller chunks, scipy supports parallelization by specifying the _workers_ environment variable:
+
+    from scipy import fft
+    with fft.set_workers(-1):
+        standardizingTrans(imgs,sigma,method,eps=0,inverse=False)
+
+When the number of workers is set to a negative integer (like above), the number of workers is set to os.cpu_count().
+
+Scipy also supports computations using another backend. For example, we can use pyFFTW as a backend like so: 
+
+     from scipy import fft
+     import pyfftw
+     with fft.set_backend(pyfftw.interfaces.scipy_fft):
+        #faster if we enable cache using pyfftw
+        pyfftw.interfaces.cache.enable()
+        # perform standardizing transform using frequency method of your choice
+        imgs_ST = standardizingTrans(imgs,sigma=(10.,20.),method='DCT',eps=0,inverse=False)
+
+## Inverse Transforms
+In principle, the frequency methods allow simple inverse ST by inverse filtering, and this can be triggered using the *inverse=True* flag. However, in the present of any noise, catastrophic noise amplification can occur, especially since the filter used in this transform is Gaussian and tends towards zero. The effect will increase with sigma. If an inverse transformation is necessary, it is recommended to use the **eps** parameter to add a small constant like 1e-5 to the filter in frequency space, i.e. the optical transfer function. The **eps** parameter is also available for (forward) standardizing transforms.
+
+

IMED.egg-info/SOURCES.txt

@@ -0,0 +1,13 @@
+LICENSE
+README.md
+pyproject.toml
+setup.cfg
+IMED/ST_all.py
+IMED/__init__.py
+IMED/benchmark.py
+IMED/spatial_ST.py
+IMED/standardizingTrans_ndim.py
+IMED.egg-info/PKG-INFO
+IMED.egg-info/SOURCES.txt
+IMED.egg-info/dependency_links.txt
+IMED.egg-info/top_level.txt

IMED.egg-info/dependency_links.txt

@@ -0,0 +1 @@
+

IMED.egg-info/top_level.txt

@@ -0,0 +1 @@
+IMED

build/lib/IMED/ST_all.py

@@ -0,0 +1,52 @@
+from IMED.standardizingTrans_ndim import ST_ndim_DCT, ST_ndim_FFT, ST_ndim_DCT_by_FFT
+from IMED.spatial_ST import ST_fullMat, ST_sepMat
+def standardizingTrans(imgs,sigma,method='dct',eps=0,inverse=False):
+    """
+    Takes n-dimensional data and returns the n-dimensional Standardized Transform.
+    Methods 'full' and 'sep' are 2D methods only.
+    
+    Parameters:
+    * imgs is a signal 
+    * sigma (float)/array-like determines the zero-mean Gaussian that defines the IMED matrix G - not G^(1/2).
+      If sigma is array-like it should contain the same number of values as the number of dimensions of imgs.
+      
+    * eps (float) is an optional small parameter to offset the Gaussian so that it is always numerically non-zero. 
+    This can allow deconvolution without significant noise amplification.
+    * method (string) is the method used to perform the standardizing transform. Choose between:
+     1. **'full':** Full Dense Matrix $z_{ST}= G^{1/2}z$ using eigenvalue decomposition
+     2. **'sep'** Separated Dense Matrices $z_{ST}= G_x^{1/2}z G_y^{1/2}$ using eigenvalue decomposition 
+     3. **'fft'**: Performs circular convolution using discrete fourier transforms of image and Gaussian 
+     without enforcing symmetric boundary conditions
+     4. **'dct_by_fft'**: Performs circular convolution using discrete fourier transforms of mirrored image and 
+     Gaussian to ensure symmetric boundary conditions and reduce edge effects from 'Gibbs-like phenomenon'
+     5. **'dct'**: Performs symmetric convolution using discrete cosine transform of image, which is identical to
+     6. the 'dct_by_fft  method, but should be more efficient
+    """
+
+    if  method == 'full':
+        if inverse==True:
+            print("No inverse method implemented")
+            return
+        imgs_ST = ST_fullMat(imgs,sigma,eps)
+
+    elif method == 'sep':
+        if inverse==True:
+            print("No inverse method implemented")
+            return
+        imgs_ST = ST_sepMat(imgs,sigma,eps)
+
+    elif method == 'fft':
+        imgs_ST = ST_ndim_FFT(imgs, sigma, eps,inverse)
+
+    elif method == 'dct_by_fft':
+        imgs_ST = ST_ndim_DCT_by_FFT(imgs, sigma, eps,inverse)
+
+    elif method == 'dct':
+        imgs_ST = ST_ndim_DCT(imgs, sigma, eps, inverse)
+
+    else:
+        print(f'Invalid method "{method}". Choosing dct.')
+        method = 'dct'
+        standardizingTrans(imgs,sigma,method,eps,inverse)    
+
+    return imgs_ST

build/lib/IMED/spatial_ST.py

@@ -0,0 +1,92 @@
+import numpy as np
+
+def ST_fullMat(imgs,sigma,eps):
+    # This is created by Niklas Heim / James Avery, and bugfixed by Jacob Felding
+    # Purposefully not M, N
+
+    #sigma check
+    if isinstance(sigma, (list, tuple, np.ndarray)):
+        # if multiple sigmas are given, there should be 2 for spatial transform
+        print('This method allows only single sigma values. Using the first one')
+        sigma = sigma[0]
+    if len(imgs.shape) == 3:
+        N, M = imgs.shape[1:]
+    else:
+        N, M = imgs.shape
+
+    X = np.arange(M,dtype=imgs.dtype)
+    Y = np.arange(N,dtype=imgs.dtype)
+
+    P = (X[None, :, None, None] - X[None, None, None, :])**2 \
+        + (Y[:, None, None, None] - Y[None, None, :, None])**2
+
+    G = 1 / (2 * np.pi * sigma**2) * np.exp(- P / (2 * sigma**2)) 
+
+    #Construct NM x NM G matrix
+    G = G.reshape((M * N, M * N))
+
+    # Use eigenvalue decomposition to construct G^1/2 used for standardizing transform
+    (w,V) = np.linalg.eigh(G)
+    # G^(1/2) has eigenvalues that are the square root of eigenvalues of G
+    w = np.maximum(np.zeros(len(w),dtype=imgs.dtype),w)
+    s = np.sqrt(w)
+
+    # Eigenvalue decomposition
+    #G_sqrt  = np.dot(V, s[:,None]*V.T)
+    VTS = s[:,None]*V.T
+
+    del G, P
+    #G_sqrt = #
+    return np.array([np.dot(V, np.dot(VTS,z.reshape(-1))) + eps for z in imgs]).reshape(-1,N,M)
+
+
+
+def ST_sepMat(imgs,sigma,eps):
+    '''Implements standardizing transform of image sequence using kronecker product matrix separation'''
+
+    # Make sigma 2-dimensional if not already
+    if isinstance(sigma, (list, tuple, np.ndarray)):
+        # Should be 2 sigmas for spatial transform
+        assert len(sigma) == 2
+    else:
+        #if sigma is a scalar, it will be used for all axes
+        sigma = np.ones(2)*sigma
+
+
+    if len(imgs.shape) == 3:
+        N, M = imgs.shape[1:]
+    else:
+        N, M = imgs.shape
+
+    # Create 1D arrays corresponding to image dimensions
+    X = np.arange(M,dtype=imgs.dtype)
+    Y = np.arange(N,dtype=imgs.dtype)  
+
+    # Create 2D arrays where P_x(i,j) denotes (x_i-x_j)**2
+    P_x = (X[None,:] - X[:,None])**2 
+    P_y = (Y[None,:] - Y[:, None])**2
+
+    # Create G decompositions such that G = np.kron(G_x,G_y)
+    # This way we store NxN, MxM matrices, not NMxNM
+    G_x =  1 / (np.sqrt((2 * np.pi ))*sigma[0]) * np.exp(- P_x / (2 * sigma[0]**2)) 
+    G_y =  1 / (np.sqrt((2 * np.pi ))*sigma[0]) * np.exp(- P_y / (2 * sigma[1]**2)) 
+
+    # Determine eigenvectors and of G using that G = np.kron(G_x,G_y)
+    # The matrices are symmetric since (x_i-x_j)**2 = (x_j-x_i)**2
+    eigvals_Gx,eigvecs_Gx = np.linalg.eigh(G_x)
+    eigvals_Gy,eigvecs_Gy = np.linalg.eigh(G_y)
+
+    # For robustness to negative eigenvalues that should not arise as
+    #G is positive definte, but may for numerical reasons
+    eigvals_Gx = np.maximum(eigvals_Gx,np.zeros(M,dtype=imgs.dtype))
+    eigvals_Gy = np.maximum(eigvals_Gy,np.zeros(N,dtype=imgs.dtype))
+
+    # For the Standardizing Transfrom G^(1/2)x (which blurs image x) 
+    # we need matrix G^(1/2) Such that G = G^(1/2)G^(1/2)
+    # Again decomposition allows G^(1/2) = np.kron(G_sqrt_x,G_sqrt_y)
+    # Below, use eigendecomposition of G to construct G_sqrt_x, G_sqrt_y
+    G_sqrt_x = np.dot(eigvecs_Gx, np.sqrt(eigvals_Gx)[:,None]*eigvecs_Gx.T) + eps
+    G_sqrt_y = np.dot(eigvecs_Gy, np.sqrt(eigvals_Gy)[:,None]*eigvecs_Gy.T) + eps
+
+    #return np.dot(G_sqrt_x,np.dot(imgs,G_sqrt_y))
+    return [np.dot(G_sqrt_x,np.dot(z,G_sqrt_y)) for z in imgs]