diff --git a/IMED.egg-info/PKG-INFO b/IMED.egg-info/PKG-INFO
new file mode 100644
index 0000000..a6dea33
--- /dev/null
+++ b/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.
+
+
diff --git a/IMED.egg-info/SOURCES.txt b/IMED.egg-info/SOURCES.txt
new file mode 100644
index 0000000..a787b4e
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/IMED.egg-info/dependency_links.txt b/IMED.egg-info/dependency_links.txt
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/IMED.egg-info/dependency_links.txt
@@ -0,0 +1 @@
+
diff --git a/IMED.egg-info/top_level.txt b/IMED.egg-info/top_level.txt
new file mode 100644
index 0000000..362ecf2
--- /dev/null
+++ b/IMED.egg-info/top_level.txt
@@ -0,0 +1 @@
+IMED
diff --git a/build/lib/IMED/ST_all.py b/build/lib/IMED/ST_all.py
new file mode 100644
index 0000000..bddd9ea
--- /dev/null
+++ b/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
diff --git a/build/lib/IMED/spatial_ST.py b/build/lib/IMED/spatial_ST.py
new file mode 100644
index 0000000..2b83393
--- /dev/null
+++ b/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]
diff --git a/build/lib/IMED/standardizingTrans_ndim.py b/build/lib/IMED/standardizingTrans_ndim.py
new file mode 100644
index 0000000..ed6a90b
--- /dev/null
+++ b/build/lib/IMED/standardizingTrans_ndim.py
@@ -0,0 +1,167 @@
+import numpy as np
+from scipy.fft import dct, idct
+#from scipy.fft import rfftn, irfftn
+#from jax.numpy.fft import rfftn, irfftn
+import jax.numpy as jnp
+from scipy.fft import rfftn, irfftn, fftfreq,rfftfreq
+
+
+NA = np.newaxis
+
+def axis_split(A,axis):
+    (Am,Ai,An) = (int(np.prod(A.shape[:axis])),  A.shape[axis], int(np.prod(A.shape[axis+1:])));
+    return (Am,Ai,An)
+
+def ST_1dim_DCT(Img,sigma,d=0,eps=0,inverse=False):
+    # Gaussian in k-space
+    k_d     = np.linspace(0,np.pi,Img.shape[d])
+    g12_dct = np.exp(- k_d**2*sigma**2/4) + eps
+
+    # Transform to k-space
+    img_dct = dct(Img, axis=d,type=1).reshape( axis_split(Img,d) )
+    
+    if inverse:
+        Img_folded_k = img_dct / g12_dct[NA,:,NA]
+    else:
+        Img_folded_k = img_dct * g12_dct[NA,:,NA]
+
+    Img_folded = idct(Img_folded_k,axis=1,type=1)
+
+    return Img_folded.reshape(Img.shape)
+    
+def ST_1dim_FFT(Img,sigma,d=0,eps=0,inverse=False,jax_backend=False):   
+    # Transform to k-space
+    #img_fft = rfftn(Img)
+    if jax_backend:
+        img_fft = jnp.fft.rfftn(Img)
+    else:   
+        img_fft = rfftn(Img)
+
+    fft_shape = img_fft.shape
+    img_fft = img_fft.reshape(axis_split(img_fft,d))
+    
+    # if last axis, need other k definition for rfft
+    if d == Img.ndim-1:
+        k_d = rfftfreq(2*fft_shape[d]-1)*2*jnp.pi
+        if Img.dtype == jnp.float32:
+            k_d = jnp.float32(k_d)
+    else:
+        k_d = fftfreq(fft_shape[d])*2*jnp.pi  
+        if Img.dtype == jnp.float32:
+            k_d = jnp.float32(k_d)
+            
+    # Gaussian in k-space
+    g12_fft = np.exp(- k_d**2*sigma**2/4) + eps
+    
+    if inverse:
+        Img_folded_k = img_fft / g12_fft[NA,:,NA]
+    else:
+        Img_folded_k = img_fft * g12_fft[NA,:,NA]
+    print()
+    if jax_backend:
+        Img_folded = jnp.fft.irfftn(Img_folded_k.reshape(fft_shape))
+    else:
+        Img_folded = irfftn(Img_folded_k.reshape(fft_shape))
+
+    return Img_folded.reshape(Img.shape)    
+    
+def ST_ndim_DCT(imgs,sigma,eps=0.,inverse=False):
+    # automatic d-dimensional standardizing transform
+    # via DCT, i.e. symmetric boundary conditions
+    # eps is an optional constant added to the OTF to reduce  
+    # noise amplification when deconvolving
+    shape = imgs.shape
+    dims     = len(shape)
+    all_dims = range(dims)
+    
+    
+    # Make sigma d-dimensional if not already
+    if isinstance(sigma, (list, tuple, np.ndarray)):
+        # if multiple sigmas are given, they must match number of axes
+        assert len(sigma) == dims
+    else:
+        #if sigma is a scalar, it will be used for all axes
+        sigma = np.ones(dims)*sigma
+        
+    # do convolution, axis by axis
+    for d in range(dims):
+        if sigma[d]==0:
+            #convolution has no effect
+            continue
+
+        if shape[d]<2:
+            #cant do convolution along this axis 
+            continue
+        
+        imgs = ST_1dim_DCT(Img=imgs,sigma=sigma[d],d=d,eps=eps,inverse=inverse)
+        
+        
+def ST_DCT_by_FFT(imgs, sigma, eps=0.,inverse=False,jax_backend=True):
+    # automatic d-dimensional standardizing transform
+    # via FFT. Uses per-axis mirroring to reduce edge discontinuities
+    # eps is an optional constant added to the OTF to reduce  
+    # noise amplification when deconvolving
+
+    orig_shape    = imgs.shape
+    dims     = len(orig_shape)
+    all_dims = range(dims)
+
+    # Make sigma d-dimensional if not already
+    if isinstance(sigma, (list, tuple, np.ndarray)):
+        # if multiple sigmas are given, they must match number of axes
+        assert len(sigma) == dims
+    else:
+        #if sigma is a scalar, it will be used for all axes
+        sigma = np.ones(dims)*sigma
+    
+    for d in range(dims):
+        if sigma[d]==0:
+            #convolution has no effect
+            continue        
+        if orig_shape[d]<3:
+            print(f'Skipping dimension d={d} with size {orig_shape[d]}')
+            #cant do convolution along this axis 
+            continue
+        
+        #mirror image along axis d
+        imgs_reverse = np.flip(imgs, axis=d)
+        # for DCT equivalence
+        imgs_reverse = imgs_reverse.take(indices=range(1,imgs_reverse.shape[d]-1),axis=d)
+        imgs         = np.concatenate((imgs,imgs_reverse),axis=d)
+        
+        imgs = ST_1dim_FFT(imgs,sigma[d],d,eps,inverse,jax_backend)
+        
+        #Cut to original shape before moving on to other axis         
+        imgs = imgs.take(indices=range(orig_shape[d]),axis=d)
+    
+    return imgs    
+ 
+def ST_ndim_FFT(imgs, sigma, eps=0.,inverse=False,jax_backend=False):
+    # automatic d-dimensional standardizing transform
+    # via FFT. Uses per-axis mirroring to reduce edge discontinuities
+    # eps is an optional constant added to the OTF to reduce  
+    # noise amplification when deconvolving
+    orig_shape    = imgs.shape
+    dims     = len(orig_shape)
+    all_dims = range(dims)
+    
+    # Make sigma d-dimensional if not already
+    if isinstance(sigma, (list, tuple, np.ndarray)):
+        # if multiple sigmas are given, they must match number of axes
+        assert len(sigma) == dims
+    else:
+        #if sigma is a scalar, it will be used for all axes
+        sigma = np.ones(dims)*sigma
+    
+    for d in range(dims):
+        if sigma[d]==0:
+            #convolution has no effect
+            continue
+            
+        if orig_shape[d]<2:
+            #cant do convolution along this axis 
+            continue
+        
+        imgs = ST_1dim_FFT(imgs,sigma[d],d,eps,inverse,jax_backend)
+    
+    return imgs
diff --git a/dist/IMED-0.0.1.8-py3-none-any.whl b/dist/IMED-0.0.1.8-py3-none-any.whl
new file mode 100644
index 0000000..1467e2a
Binary files /dev/null and b/dist/IMED-0.0.1.8-py3-none-any.whl differ
diff --git a/dist/IMED-0.0.1.8.tar.gz b/dist/IMED-0.0.1.8.tar.gz
new file mode 100644
index 0000000..c0c9f98
Binary files /dev/null and b/dist/IMED-0.0.1.8.tar.gz differ
diff --git a/setup.cfg b/setup.cfg
index 6dd98f3..859f9c4 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,7 +1,7 @@
 [metadata]
 # replace with your username:
 name = IMED
-version = 0.0.1.7-a
+version = 0.0.1.8
 author = Jacob Felding
 author_email = jfelding@gmail.com
 description = Image Euclidean Distance and Standardizing Transforms