C lib (#458)

C library implemented with optimised axpy fucntions and gradient operator in c

23 files changed, 1695 insertions, 38 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100755
index 0000000..2df50b5
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,7 @@
+cmake_minimum_required(VERSION 3.4)
+
+project (cil LANGUAGES C CXX)
+
+set(CIL_VERSION $ENV{CIL_VERSION})
+add_subdirectory(src/Core)
+add_subdirectory(Wrappers/Python)
diff --git a/Wrappers/Python/CMakeLists.txt b/Wrappers/Python/CMakeLists.txt
new file mode 100644
index 0000000..f325265
--- /dev/null
+++ b/Wrappers/Python/CMakeLists.txt
@@ -0,0 +1,88 @@
+option (BUILD_PYTHON_WRAPPER "Build Python Wrapper" ON)
+
+if (BUILD_PYTHON_WRAPPER)
+  find_package(PythonInterp REQUIRED)
+  
+  set(PYTHON_DEST_DIR "" CACHE PATH "Directory of the Python wrappers")
+  if (PYTHON_DEST_DIR)
+  set(PYTHON_DEST "${PYTHON_DEST_DIR}")
+  else() 
+  set(PYTHON_DEST "${CMAKE_INSTALL_PREFIX}/python")
+  endif()
+  message(STATUS "Python wrappers will be installed in " ${PYTHON_DEST})
+
+  message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}")
+
+  set(CMAKE_BUILD_TYPE "Release")
+  
+  find_package(PythonLibs)
+  if (PYTHONINTERP_FOUND)
+    message(STATUS "Found PYTHON_EXECUTABLE=${PYTHON_EXECUTABLE}")
+    message(STATUS "Python version ${PYTHON_VERSION_STRING}")
+  endif()
+  if (PYTHONLIBS_FOUND)
+    message(STATUS "Found PYTHON_INCLUDE_DIRS=${PYTHON_INCLUDE_DIRS}")
+    message(STATUS "Found PYTHON_LIBRARIES=${PYTHON_LIBRARIES}")
+  endif()
+
+  if (PYTHONINTERP_FOUND)
+      message("Python found " ${PYTHON_EXECUTABLE})
+      set(SETUP_PY_IN "${CMAKE_CURRENT_SOURCE_DIR}/setup-cil.py.in")
+      set(SETUP_PY    "${CMAKE_CURRENT_BINARY_DIR}/setup.py")
+      #set(DEPS        "${CMAKE_CURRENT_SOURCE_DIR}/module/__init__.py")
+      set (DEPS       "${CMAKE_BINARY_DIR}/")
+      set(OUTPUT      "${CMAKE_CURRENT_BINARY_DIR}/build/timestamp")
+
+      #configure_file(${SETUP_PY_IN} ${SETUP_PY})
+
+      message("Core binary dir " ${CMAKE_BINARY_DIR}/Core/${CMAKE_BUILD_TYPE})
+      
+      if (CONDA_BUILD)
+        add_custom_command(OUTPUT ${OUTPUT}
+                        #COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}
+                        COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}
+                        COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION}
+                                                        #PREFIX=${CMAKE_SOURCE_DIR}/src/Core
+                                                        #LIBRARY_INC=${CMAKE_SOURCE_DIR}/src/Core
+                                                        #LIBRARY_LIB=${CMAKE_BINARY_DIR}/src/Core
+                                                        ${PYTHON_EXECUTABLE} ${SETUP_PY} -vv install
+                                                        #echo "EDO"
+                        WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+                        COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
+                        DEPENDS cilacc)
+
+      else()
+        if (WIN32)
+          add_custom_command(OUTPUT ${OUTPUT}
+                        COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}
+                        COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION}
+                                                        PREFIX=${CMAKE_SOURCE_DIR}/src/
+                                                        LIBRARY_INC=${CMAKE_SOURCE_DIR}/src/include
+                                                        LIBRARY_LIB=${CMAKE_BINARY_DIR}/${CMAKE_BUILD_TYPE}
+                                                        ${PYTHON_EXECUTABLE} ${SETUP_PY} build_ext --inplace
+                        COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
+                        DEPENDS cilacc)
+        else()
+          add_custom_command(OUTPUT ${OUTPUT}
+                        COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}
+                        COMMAND ${CMAKE_COMMAND} -E env PREFIX=${CMAKE_SOURCE_DIR}/src/
+                                                        LIBRARY_INC=${CMAKE_SOURCE_DIR}/src/include
+                                                        LIBRARY_LIB=${CMAKE_BINARY_DIR}/
+                                                        ${PYTHON_EXECUTABLE} ${SETUP_PY} build_ext --inplace
+                        COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
+                        DEPENDS cilacc
+                        )
+        endif()
+        #set (PYTHON_DEST ${CMAKE_INSTALL_PREFIX}/python/)
+        install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/
+                DESTINATION ${PYTHON_DEST})
+        #file(TOUCH ${PYTHON_DEST}/edo/__init__.py)
+        
+      endif()
+      
+      
+      add_custom_target(PythonWrapper ALL DEPENDS ${OUTPUT})
+      #install(CODE "execute_process(COMMAND ${PYTHON} ${SETUP_PY} install)")
+  endif()
+
+endif()
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py
index 01ce7ef..1ab1d1e 100755
--- a/Wrappers/Python/ccpi/framework/framework.py
+++ b/Wrappers/Python/ccpi/framework/framework.py
@@ -26,6 +26,26 @@ from datetime import timedelta, datetime
 import warnings
 from functools import reduce
 from numbers import Number
+import ctypes, platform
+
+# dll = os.path.abspath(os.path.join( 
+#          os.path.abspath(os.path.dirname(__file__)),
+#          'libfdiff.dll')
+# )
+
+# check for the extension
+if platform.system() == 'Linux':
+    dll = 'libcilacc.so'
+elif platform.system() == 'Windows':
+    dll = 'cilacc.dll'
+elif platform.system() == 'Darwin':
+    dll = 'libcilacc.dylib'
+else:
+    raise ValueError('Not supported platform, ', platform.system())
+
+#print ("dll location", dll)
+cilacc = ctypes.cdll.LoadLibrary(dll)
+
 
 def find_key(dic, val):
     """return the key of dictionary dic given the value"""
@@ -802,7 +822,69 @@ class DataContainer(object):
     def minimum(self,x2, out=None, *args, **kwargs):
         return self.pixel_wise_binary(numpy.minimum, x2=x2, out=out, *args, **kwargs)
 
-    
+    @staticmethod
+    def axpby(a,x,b,y,out,dtype=numpy.float32):
+        '''performs axpby with cilacc C library
+        
+        Does the operation .. math:: a*x+b*y and stores the result in out
+
+        :param a: scalar
+        :param x: DataContainer
+        :param b: scalar
+        :param y: DataContainer
+        :param out: DataContainer to store the result
+        :param dtype: optional, data type of the DataContainers
+        '''
+
+        c_float_p = ctypes.POINTER(ctypes.c_float)
+        c_double_p = ctypes.POINTER(ctypes.c_double)
+        # get the reference to the data
+        ndx = x.as_array()
+        ndy = y.as_array()
+        ndout = out.as_array()
+
+        if ndx.dtype != dtype:
+            ndx = ndx.astype(dtype)
+        if ndy.dtype != dtype:
+            ndy = ndy.astype(dtype)
+        
+        if dtype == numpy.float32:
+            x_p = ndx.ctypes.data_as(c_float_p)
+            y_p = ndy.ctypes.data_as(c_float_p)
+            out_p = ndout.ctypes.data_as(c_float_p)
+            f = cilacc.saxpby
+
+        elif dtype == numpy.float64:
+            ndx = ndx.astype(numpy.float64)
+            b = b.astype(numpy.float64)
+            x_p = ndx.ctypes.data_as(c_double_p)
+            y_p = ndy.ctypes.data_as(c_double_p)
+            out_p = ndout.ctypes.data_as(c_double_p)
+            f = cilacc.daxpby
+        else:
+            raise TypeError('Unsupported type {}. Expecting numpy.float32 or numpy.float64'.format(dtype))
+
+        #out = numpy.empty_like(a)
+
+        
+        # int psaxpby(float * x, float * y, float * out, float a, float b, long size)
+        cilacc.saxpby.argtypes = [ctypes.POINTER(ctypes.c_float),  # pointer to the first array 
+                                  ctypes.POINTER(ctypes.c_float),  # pointer to the second array 
+                                  ctypes.POINTER(ctypes.c_float),  # pointer to the third array 
+                                  ctypes.c_float,                  # type of A (float)
+                                  ctypes.c_float,                  # type of B (float)
+                                  ctypes.c_long]                   # type of size of first array 
+        cilacc.daxpby.argtypes = [ctypes.POINTER(ctypes.c_double), # pointer to the first array 
+                                  ctypes.POINTER(ctypes.c_double), # pointer to the second array 
+                                  ctypes.POINTER(ctypes.c_double), # pointer to the third array 
+                                  ctypes.c_double,                 # type of A (c_double)
+                                  ctypes.c_double,                 # type of B (c_double)
+                                  ctypes.c_long]                   # type of size of first array 
+
+        if f(x_p, y_p, out_p, a, b, ndx.size) != 0:
+            raise RuntimeError('axpby execution failed')
+        
+
     ## unary operations
     def pixel_wise_unary(self, pwop, *args,  **kwargs):
         out = kwargs.get('out', None)
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index 3c32a93..8e07802 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -21,15 +21,37 @@ from __future__ import print_function
 from __future__ import unicode_literals
 
 from ccpi.optimisation.operators import Operator, LinearOperator, ScaledOperator
+from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff
 from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer
 import numpy 
-from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff
+import warnings
 
-#%%
+NEUMANN = 'Neumann'
+PERIODIC = 'Periodic'
+C = 'c'
+NUMPY = 'numpy'
+CORRELATION_SPACE = "Space"
+CORRELATION_SPACECHANNEL = "SpaceChannels"
 
 class Gradient(LinearOperator):
+
+    """This is a class to compute the first-order forward/backward differences on ImageData
     
+    :param gm_domain: Set up the domain of the function
+    :type gm_domain: `ImageGeometry`
+    :param bnd_cond: Set the boundary conditions to use 'Neumann' or 'Periodic', defaults to 'Neumann'
+    :type bnd_cond: str, optional    
     
+    :param \**kwargs:
+        See below
+
+    :Keyword Arguments:
+        * *correlation* (``str``) --
+          'Space' or 'SpaceChannels', defaults to 'Space'
+        * *backend* (``str``) --
+          'c' or 'numpy', defaults to 'c' if correlation is 'SpaceChannels' or channels = 1
+    """
+
     r'''Gradient Operator: .. math:: \nabla : X -> Y           
             
             Computes first-order forward/backward differences 
@@ -39,29 +61,82 @@ class Gradient(LinearOperator):
                 Example (2D): u\in X, \nabla(u) = [\partial_{y} u, \partial_{x} u]
                               u^{*}\in Y, \nabla^{*}(u^{*}) = \partial_{y} v1 + \partial_{x} v2
         
-                Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
-                Grad_order = ['channels', 'direction_y', 'direction_x']
-                Grad_order = ['direction_z', 'direction_y', 'direction_x']
-                Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
-                
-
     '''
+
+    #kept here for backwards compatability
+    CORRELATION_SPACE = CORRELATION_SPACE
+    CORRELATION_SPACECHANNEL = CORRELATION_SPACECHANNEL
+
+    def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs):
+        """Constructor method
+        """        
+        super(Gradient, self).__init__() 
+
+        backend = kwargs.get('backend',C)
+
+        correlation = kwargs.get('correlation',CORRELATION_SPACE)
+
+        if correlation == CORRELATION_SPACE and gm_domain.channels > 1:
+            #numpy implementation only for now
+            backend = NUMPY
+            warnings.warn("Warning: correlation='Space' on multi-channel dataset will use `numpy` backend")
+
+        if backend == NUMPY:
+            self.operator = Gradient_numpy(gm_domain, bnd_cond=bnd_cond, **kwargs)
+        else:
+            self.operator = Gradient_C(gm_domain, bnd_cond=bnd_cond)
+
+
+    def direct(self, x, out=None):
+        """Computes the first-order forward differences
+
+        :param x: Image data
+        :type x: `ImageData`
+        :param out: pre-allocated output memory to store result
+        :type out: `BlockDataContainer`, optional        
+        :return: result data if not passed as parameter
+        :rtype: `BlockDataContainer`
+        """        
+        return self.operator.direct(x, out=out)
+        
+        
+    def adjoint(self, x, out=None):
+        """Computes the first-order backward differences
+
+        :param x: Gradient images for each dimension in ImageGeometry domain
+        :type x: `BlockDataContainer`
+        :param out: pre-allocated output memory to store result
+        :type out: `ImageData`, optional      
+        :return: result data if not passed as parameter
+        :rtype: `ImageData`
+        """            
+        return self.operator.adjoint(x, out=out)
+
+    def domain_geometry(self):
+        '''Returns domain_geometry of Gradient'''
+        
+        return self.operator.gm_domain
     
-                    
-    CORRELATION_SPACE = "Space"
-    CORRELATION_SPACECHANNEL = "SpaceChannels"
+    def range_geometry(self):
+        
+        '''Returns range_geometry of Gradient'''
+        
+        return self.operator.gm_range
 
+class Gradient_numpy(LinearOperator):
+    
     def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs):
         
-        super(Gradient, self).__init__() 
+        super(Gradient_numpy, self).__init__() 
                 
         self.gm_domain = gm_domain # Domain of Grad Operator
         
-        self.correlation = kwargs.get('correlation',Gradient.CORRELATION_SPACE)
+        self.correlation = kwargs.get('correlation',CORRELATION_SPACE)
         
-        if self.correlation==Gradient.CORRELATION_SPACE:
-            if self.gm_domain.channels>1:
+        if self.correlation==CORRELATION_SPACE:
+            if self.gm_domain.channels > 1:
                 self.gm_range = BlockGeometry(*[self.gm_domain for _ in range(self.gm_domain.length-1)] )
+
                 if self.gm_domain.length == 4:
                     # 3D + Channel
                     # expected Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
@@ -70,6 +145,7 @@ class Gradient(LinearOperator):
                     # 2D + Channel
                     # expected Grad_order = ['channels', 'direction_y', 'direction_x']
                     expected_order = [ImageGeometry.CHANNEL, ImageGeometry.HORIZONTAL_Y, ImageGeometry.HORIZONTAL_X]
+
                 order = self.gm_domain.get_order_by_label(self.gm_domain.dimension_labels, expected_order)
                 self.ind = order[1:]
                 #self.ind = numpy.arange(1,self.gm_domain.length)
@@ -83,10 +159,12 @@ class Gradient(LinearOperator):
                 else:
                     # 2D
                     expected_order = [ImageGeometry.HORIZONTAL_Y, ImageGeometry.HORIZONTAL_X]    
+
                 self.ind = self.gm_domain.get_order_by_label(self.gm_domain.dimension_labels, expected_order)
                 # self.ind = numpy.arange(self.gm_domain.length)
-        elif self.correlation==Gradient.CORRELATION_SPACECHANNEL:
-            if self.gm_domain.channels>1:
+                
+        elif self.correlation==CORRELATION_SPACECHANNEL:
+            if self.gm_domain.channels > 1:
                 self.gm_range = BlockGeometry(*[self.gm_domain for _ in range(self.gm_domain.length)])
                 self.ind = range(self.gm_domain.length)
             else:
@@ -189,6 +267,146 @@ class Gradient(LinearOperator):
             res.append(spMat.sum_abs_col())
         return BlockDataContainer(*res)
    
+import ctypes, platform
+
+# check for the extension
+if platform.system() == 'Linux':
+    dll = 'libcilacc.so'
+elif platform.system() == 'Windows':
+    dll = 'cilacc.dll'
+elif platform.system() == 'Darwin':
+    dll = 'libcilacc.dylib'
+else:
+    raise ValueError('Not supported platform, ', platform.system())
+
+#print ("dll location", dll)
+cilacc = ctypes.cdll.LoadLibrary(dll)
+
+#FD = ctypes.CDLL(dll)
+
+c_float_p = ctypes.POINTER(ctypes.c_float)
+
+cilacc.openMPtest.restypes = ctypes.c_int32
+
+cilacc.fdiff4D.argtypes = [ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.c_long,
+                       ctypes.c_long,
+                       ctypes.c_long,
+                       ctypes.c_long,
+                       ctypes.c_int32,
+                       ctypes.c_int32]
+
+cilacc.fdiff3D.argtypes = [ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.c_long,
+                       ctypes.c_long,
+                       ctypes.c_long,
+                       ctypes.c_int32,
+                       ctypes.c_int32]
+
+cilacc.fdiff2D.argtypes = [ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.c_long,
+                       ctypes.c_long,
+                       ctypes.c_int32,
+                       ctypes.c_int32]
+
+
+class Gradient_C(LinearOperator):
+    
+    '''Finite Difference Operator:
+            
+            Computes first-order forward/backward differences 
+                     on 2D, 3D, 4D ImageData
+                     under Neumann/Periodic boundary conditions'''
+
+    def __init__(self, gm_domain, gm_range=None, bnd_cond = NEUMANN):
+
+        super(Gradient_C, self).__init__() 
+
+        self.gm_domain = gm_domain
+        self.gm_range = gm_range
+        
+        #default is 'Neumann'
+        self.bnd_cond = 0
+        
+        if bnd_cond == PERIODIC:
+            self.bnd_cond = 1
+        
+        # Domain Geometry = Range Geometry if not stated
+        if self.gm_range is None:
+            self.gm_range = BlockGeometry(*[gm_domain for _ in range(len(gm_domain.shape))])
+        
+
+        if len(gm_domain.shape) == 4:
+            self.fd = cilacc.fdiff4D
+        elif len(gm_domain.shape) == 3:
+            self.fd = cilacc.fdiff3D
+        elif len(gm_domain.shape) == 2:
+            self.fd = cilacc.fdiff2D
+        else:
+            raise ValueError('Number of dimensions not supported, expected 2, 3 or 4, got {}'.format(len(gm_domain.shape)))
+        
+                        
+    @staticmethod 
+    def datacontainer_as_c_pointer(x):
+        ndx = x.as_array()
+        return ndx, ndx.ctypes.data_as(c_float_p)
+        
+    def direct(self, x, out=None):
+        ndx , x_p = Gradient_C.datacontainer_as_c_pointer(x)
+        
+        return_val = False
+        if out is None:
+            out = self.gm_range.allocate(None)
+            return_val = True
+
+        arg1 = [Gradient_C.datacontainer_as_c_pointer(out.get_item(i))[1] for i in range(self.gm_range.shape[0])]
+        arg2 = [el for el in x.shape]
+        args = arg1 + arg2 + [self.bnd_cond, 1]
+        self.fd(x_p, *args)
+        
+        if return_val is True:
+            return out
+
+
+    def adjoint(self, x, out=None):
+
+        return_val = False
+        if out is None:
+            out = self.gm_domain.allocate(None)
+            return_val = True
+
+        ndout, out_p = Gradient_C.datacontainer_as_c_pointer(out)
+
+        arg1 = [Gradient_C.datacontainer_as_c_pointer(x.get_item(i))[1] for i in range(self.gm_range.shape[0])]
+        arg2 = [el for el in out.shape]
+        args = arg1 + arg2 + [self.bnd_cond, 0]
+
+        self.fd(out_p, *args)
+
+        if return_val is True:
+            return out
+
+    def domain_geometry(self):
+        
+        '''Returns domain_geometry of Gradient'''
+        
+        return self.gm_domain
+    
+    def range_geometry(self):
+        
+        '''Returns range_geometry of Gradient'''
+        
+        return self.gm_range
+
        
 if __name__ == '__main__':
     
@@ -302,8 +520,5 @@ if __name__ == '__main__':
     G.direct(u, out = res)          
     G.adjoint(w, out = res1)
           
-    print( (res*w).sum(),  (u*res1).sum() )     
-    
+    print( (res*w).sum(),  (u*res1).sum() )
     
-    
-
diff --git a/Wrappers/Python/conda-recipe/bld.bat b/Wrappers/Python/recipe/bld.bat
index 97a4e62..97a4e62 100644
--- a/Wrappers/Python/conda-recipe/bld.bat
+++ b/Wrappers/Python/recipe/bld.bat
diff --git a/Wrappers/Python/conda-recipe/build.sh b/Wrappers/Python/recipe/build.sh
index 43e85d5..43e85d5 100644
--- a/Wrappers/Python/conda-recipe/build.sh
+++ b/Wrappers/Python/recipe/build.sh
diff --git a/Wrappers/Python/conda-recipe/conda_build_config.yaml b/Wrappers/Python/recipe/conda_build_config.yaml
index e5e168b..e5e168b 100644
--- a/Wrappers/Python/conda-recipe/conda_build_config.yaml
+++ b/Wrappers/Python/recipe/conda_build_config.yaml
diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/recipe/meta.yaml
index 9d03220..9d03220 100644
--- a/Wrappers/Python/conda-recipe/meta.yaml
+++ b/Wrappers/Python/recipe/meta.yaml
diff --git a/Wrappers/Python/setup-cil.py.in b/Wrappers/Python/setup-cil.py.in
new file mode 100644
index 0000000..d1cf7fc
--- /dev/null
+++ b/Wrappers/Python/setup-cil.py.in
@@ -0,0 +1,70 @@
+# -*- coding: utf-8 -*-
+#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
+
+#   Copyright 2017 UKRI-STFC
+#   Copyright 2017 University of Manchester
+
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+
+#   http://www.apache.org/licenses/LICENSE-2.0
+
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+from distutils.core import setup
+import os
+import sys
+
+
+cil_version=os.environ['CIL_VERSION']
+if  cil_version == '':
+    print("Please set the environmental variable CIL_VERSION")
+    sys.exit(1)
+
+setup(
+    name="ccpi-framework",
+    version=cil_version,
+    packages=['ccpi' , 'ccpi.io',
+              'ccpi.framework', 'ccpi.optimisation', 
+              'ccpi.optimisation.operators',
+              'ccpi.optimisation.algorithms',
+              'ccpi.optimisation.functions',
+              'ccpi.processors',
+              'ccpi.utilities', 'ccpi.utilities.jupyter',
+              'ccpi.contrib','ccpi.contrib.optimisation',
+              'ccpi.contrib.optimisation.algorithms'],
+    data_files = [('share/ccpi', ['data/boat.tiff',
+                                  'data/peppers.tiff',
+                                  'data/camera.png',
+                                  'data/resolution_chart.tiff',
+                                  'data/shapes.png',
+                                  'data/24737_fd_normalised.nxs'])],
+
+    # Project uses reStructuredText, so ensure that the docutils get
+    # installed or upgraded on the target machine
+    #install_requires=['docutils>=0.3'],
+
+#    package_data={
+#        # If any package contains *.txt or *.rst files, include them:
+#        '': ['*.txt', '*.rst'],
+#        # And include any *.msg files found in the 'hello' package, too:
+#        'hello': ['*.msg'],
+#    },
+    # zip_safe = False,
+
+    # metadata for upload to PyPI
+    author="CCPi developers",
+    maintainer="Edoardo Pasca",
+    maintainer_email="edoardo.pasca@stfc.ac.uk",
+    description='CCPi Core Imaging Library - Python Framework Module',
+    license="Apache v2.0",
+    keywords="Python Framework",
+    url="http://www.ccpi.ac.uk/CIL",   # project home page, if any
+
+    # could also include long_description, download_url, classifiers, etc.
+)
diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 9b9a8dd..d2d5b05 100755
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -732,6 +732,19 @@ class TestDataContainer(unittest.TestCase):
         u.multiply(2, out=u)
         c = b * 2
         numpy.testing.assert_array_equal(u.as_array(), c)
+
+    def test_axpby(self):
+        print ("test axpby")
+        ig = ImageGeometry(10,10)                                               
+        d1 = ig.allocate(1)                                                     
+        d2 = ig.allocate(2)                                                     
+        out = ig.allocate(None)                                                 
+        # equals to 2 * [1] + 1 * [2] = [4]
+        DataContainer.axpby(2,d1,1,d2,out)
+        res = numpy.ones_like(d1.as_array()) * 4.
+        numpy.testing.assert_array_equal(res, out.as_array())
+
+
         
 
 
diff --git a/Wrappers/Python/test/test_Gradient.py b/Wrappers/Python/test/test_Gradient.py
index 5dc8137..5aeede0 100755
--- a/Wrappers/Python/test/test_Gradient.py
+++ b/Wrappers/Python/test/test_Gradient.py
@@ -23,6 +23,7 @@ from ccpi.framework import BlockDataContainer
 import functools
 
 from ccpi.optimisation.operators import Gradient, Identity, BlockOperator
+from ccpi.optimisation.operators import LinearOperator
 
 class TestGradient(unittest.TestCase):
     def test_Gradient(self): 
@@ -36,17 +37,17 @@ class TestGradient(unittest.TestCase):
         ig3 = ImageGeometry(voxel_num_x = M, voxel_num_y = N, channels = channels) 
         ig4 = ImageGeometry(voxel_num_x = M, voxel_num_y = N, channels = channels, voxel_num_z= K) 
         
-        G1 = Gradient(ig1, correlation = 'Space') 
+        G1 = Gradient(ig1, correlation = 'Space', backend='numpy') 
         print(G1.range_geometry().shape, '2D no channels')
             
-        G4 = Gradient(ig3, correlation = 'SpaceChannels')
-        print(G4.range_geometry().shape, '2D with channels corr')
-        G5 = Gradient(ig3, correlation = 'Space')
+        G4 = Gradient(ig3, correlation = 'SpaceChannels', backend='numpy')
+        print(G4.range_geometry().shape, '2D with channels corr')  
+        G5 = Gradient(ig3, correlation = 'Space', backend='numpy')
         print(G5.range_geometry().shape, '2D with channels no corr')
         
-        G6 = Gradient(ig4, correlation = 'Space')
+        G6 = Gradient(ig4, correlation = 'Space', backend='numpy')
         print(G6.range_geometry().shape, '3D with channels no corr')
-        G7 = Gradient(ig4, correlation = 'SpaceChannels')
+        G7 = Gradient(ig4, correlation = 'SpaceChannels', backend='numpy')
         print(G7.range_geometry().shape, '3D with channels with corr')
         
         
@@ -77,8 +78,8 @@ class TestGradient(unittest.TestCase):
         #check direct/adjoint for space/channels correlation
         
         ig_channel = ImageGeometry(voxel_num_x = 2, voxel_num_y = 3, channels = 2)
-        G_no_channel = Gradient(ig_channel, correlation = 'Space')
-        G_channel = Gradient(ig_channel, correlation = 'SpaceChannels')
+        G_no_channel = Gradient(ig_channel, correlation = 'Space', backend='numpy')
+        G_channel = Gradient(ig_channel, correlation = 'SpaceChannels', backend='numpy')
         
         u3 = ig_channel.allocate('random_int')
         res_no_channel = G_no_channel.direct(u3)
@@ -95,7 +96,7 @@ class TestGradient(unittest.TestCase):
         
         ig2D = ImageGeometry(voxel_num_x = 2, voxel_num_y = 3)
         u4 = ig2D.allocate('random_int')
-        G2D = Gradient(ig2D)
+        G2D = Gradient(ig2D, backend='numpy')
         res = G2D.direct(u4)  
         print(res[0].as_array())
         print(res[1].as_array())
@@ -105,12 +106,133 @@ class TestGradient(unittest.TestCase):
         arr = ig.allocate('random_int' )
         
         # check direct of Gradient and sparse matrix
-        G = Gradient(ig)
+        G = Gradient(ig, backend='numpy')
         norm1 = G.norm(iterations=300)
         print ("should be sqrt(8) {} {}".format(numpy.sqrt(8), norm1))
         numpy.testing.assert_almost_equal(norm1, numpy.sqrt(8), decimal=1)
         ig4 = ImageGeometry(M,N, channels=3)
-        G4 = Gradient(ig4, correlation=Gradient.CORRELATION_SPACECHANNEL)
+        G4 = Gradient(ig4, correlation="SpaceChannels", backend='numpy')
         norm4 = G4.norm(iterations=300)
-        print ("should be sqrt(12) {} {}".format(numpy.sqrt(12), norm4))
+        print("should be sqrt(12) {} {}".format(numpy.sqrt(12), norm4))
         self.assertTrue((norm4 - numpy.sqrt(12))/norm4 < 0.2)
+
+    def test_GradientOperator_4D(self):
+
+        nc, nz, ny, nx = 3, 4, 5, 6
+        size = nc * nz * ny * nx
+        dim = [nc, nz, ny, nx]
+
+        ig = ImageGeometry(voxel_num_x=nx, voxel_num_y=ny, voxel_num_z=nz, channels=nc)
+
+        arr = numpy.arange(size).reshape(dim).astype(numpy.float32)**2
+
+        data = ig.allocate()
+        data.fill(arr)
+
+        #neumann
+        grad_py = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='numpy')
+        gold_direct = grad_py.direct(data)
+        gold_adjoint = grad_py.adjoint(gold_direct)
+
+        grad_c = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='c')
+        out_direct = grad_c.direct(data)
+        out_adjoint = grad_c.adjoint(out_direct)
+
+        #print("GradientOperator, 4D, bnd_cond='Neumann', direct")
+        numpy.testing.assert_array_equal(out_direct.get_item(0).as_array(), gold_direct.get_item(0).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(1).as_array(), gold_direct.get_item(1).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(2).as_array(), gold_direct.get_item(2).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(3).as_array(), gold_direct.get_item(3).as_array())
+
+        #print("GradientOperator, 4D, bnd_cond='Neumann', adjoint")
+        numpy.testing.assert_array_equal(out_adjoint.as_array(), gold_adjoint.as_array())
+
+        #periodic
+        grad_py = Gradient(ig, bnd_cond='Periodic', correlation='SpaceChannels', backend='numpy')
+        gold_direct = grad_py.direct(data)
+        gold_adjoint = grad_py.adjoint(gold_direct)
+
+        grad_c = Gradient(ig, bnd_cond='Periodic', correlation='SpaceChannels', backend='c')
+        out_direct = grad_c.direct(data)
+        out_adjoint = grad_c.adjoint(out_direct)
+
+        #print("Gradient, 4D, bnd_cond='Periodic', direct")
+        numpy.testing.assert_array_equal(out_direct.get_item(0).as_array(), gold_direct.get_item(0).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(1).as_array(), gold_direct.get_item(1).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(2).as_array(), gold_direct.get_item(2).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(3).as_array(), gold_direct.get_item(3).as_array())
+
+        #print("Gradient, 4D, bnd_cond='Periodic', adjoint")
+        numpy.testing.assert_array_equal(out_adjoint.as_array(), gold_adjoint.as_array())
+
+    def test_Gradient_4D_allocate(self):
+
+        nc, nz, ny, nx = 3, 4, 5, 6
+        size = nc * nz * ny * nx
+        dim = [nc, nz, ny, nx]
+
+        ig = ImageGeometry(voxel_num_x=nx, voxel_num_y=ny, voxel_num_z=nz, channels=nc)
+
+
+        arr = numpy.arange(size).reshape(dim).astype(numpy.float32)**2
+
+        data = ig.allocate()
+        data.fill(arr)
+
+        #numpy
+        grad1 = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='numpy')
+        gold_direct = grad1.direct(data)
+        gold_adjoint = grad1.adjoint(gold_direct)
+
+        grad2 = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='numpy')
+        out_direct = grad2.range_geometry().allocate()
+        out_adjoint = grad2.domain_geometry().allocate()
+        grad2.direct(data, out=out_direct)
+        grad2.adjoint(out_direct, out=out_adjoint)
+
+        #print("GradientOperator, 4D, bnd_cond='Neumann', direct")
+        numpy.testing.assert_array_equal(out_direct.get_item(0).as_array(), gold_direct.get_item(0).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(1).as_array(), gold_direct.get_item(1).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(2).as_array(), gold_direct.get_item(2).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(3).as_array(), gold_direct.get_item(3).as_array())
+
+        #print("GradientOperator, 4D, bnd_cond='Neumann', adjoint")
+        numpy.testing.assert_array_equal(out_adjoint.as_array(), gold_adjoint.as_array())
+
+        #c
+        grad1 = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='c')
+        gold_direct = grad1.direct(data)
+        gold_adjoint = grad1.adjoint(gold_direct)
+
+        grad2 = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='c')
+        out_direct = grad2.range_geometry().allocate()
+        out_adjoint = grad2.domain_geometry().allocate()
+        grad2.direct(data, out=out_direct)
+        grad2.adjoint(out_direct, out=out_adjoint)
+
+        #print("GradientOperator, 4D, bnd_cond='Neumann', direct")
+        numpy.testing.assert_array_equal(out_direct.get_item(0).as_array(), gold_direct.get_item(0).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(1).as_array(), gold_direct.get_item(1).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(2).as_array(), gold_direct.get_item(2).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(3).as_array(), gold_direct.get_item(3).as_array())
+
+        #print("GradientOperator, 4D, bnd_cond='Neumann', adjoint")
+        numpy.testing.assert_array_equal(out_adjoint.as_array(), gold_adjoint.as_array())
+
+    def test_Gradient_linearity(self):
+
+        nc, nz, ny, nx = 3, 4, 5, 6
+        ig = ImageGeometry(voxel_num_x=nx, voxel_num_y=ny, voxel_num_z=nz, channels=nc)
+
+        grad = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='c')
+        self.assertTrue(LinearOperator.dot_test(grad))
+
+        grad = Gradient(ig, bnd_cond='Periodic', correlation='SpaceChannels', backend='c')
+        self.assertTrue(LinearOperator.dot_test(grad))
+
+        grad = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='numpy')
+        self.assertTrue(LinearOperator.dot_test(grad))
+
+        grad = Gradient(ig, bnd_cond='Periodic', correlation='SpaceChannels', backend='numpy')
+        self.assertTrue(LinearOperator.dot_test(grad))
+        
\ No newline at end of file
diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py
index 3372b9b..96c48da 100644
--- a/Wrappers/Python/test/test_Operator.py
+++ b/Wrappers/Python/test/test_Operator.py
@@ -494,12 +494,13 @@ class TestBlockOperator(unittest.TestCase):
             print("Z1", Z1[0][1].as_array())
             print("RES1", RES1[0][1].as_array())
     def test_timedifference(self):
+
         print ("test_timedifference")
         M, N ,W = 100, 512, 512
         ig = ImageGeometry(M, N, W)
         arr = ig.allocate('random_int')  
         
-        G = Gradient(ig)
+        G = Gradient(ig, backend='numpy')
         Id = Identity(ig)
         
         B = BlockOperator(G, Id)
diff --git a/build/jenkins-build.sh b/build/jenkins-build.sh
deleted file mode 100644
index 009d43d..0000000
--- a/build/jenkins-build.sh
+++ /dev/null
@@ -1,3 +0,0 @@
-#!/usr/bin/env bash
-export CCPI_BUILD_ARGS="-c conda-forge -c ccpi"
-bash <(curl -L https://raw.githubusercontent.com/vais-ral/CCPi-VirtualMachine/master/scripts/jenkins-build.sh)
diff --git a/recipe/bld.bat b/recipe/bld.bat
new file mode 100644
index 0000000..af5ca40
--- /dev/null
+++ b/recipe/bld.bat
@@ -0,0 +1,15 @@
+
+IF NOT DEFINED CIL_VERSION (
+ECHO CIL_VERSION Not Defined.
+exit 1
+)
+
+ROBOCOPY /E "%RECIPE_DIR%\.." "%SRC_DIR%" /XD .git /XD Wrappers\Python\build
+
+mkdir "%SRC_DIR%\build_framework"
+
+cd "%SRC_DIR%\build_framework"
+cmake -G "NMake Makefiles" %RECIPE_DIR%\..\ -DCONDA_BUILD=ON -DCMAKE_BUILD_TYPE="Release" -DLIBRARY_LIB=%CONDA_PREFIX%\lib -DLIBRARY_INC=%CONDA_PREFIX% -DCMAKE_INSTALL_PREFIX=%PREFIX%
+
+nmake install
+if errorlevel 1 exit 1
diff --git a/recipe/build.sh b/recipe/build.sh
new file mode 100644
index 0000000..22cfae8
--- /dev/null
+++ b/recipe/build.sh
@@ -0,0 +1,24 @@
+if [ -z "$CIL_VERSION" ]; then
+    echo "Need to set CIL_VERSION"
+    exit 1
+fi
+# mkdir ${SRC_DIR}/ccpi
+mkdir -p ${SRC_DIR}/ccpi/Wrappers/Python
+cp -r "${RECIPE_DIR}/../Wrappers/Python/test" ${SRC_DIR}/ccpi/Wrappers/Python
+
+# cd ${SRC_DIR}/ccpi/Wrappers/Python
+# $PYTHON setup.py install
+
+
+mkdir ${SRC_DIR}/build_framework
+#cp -r "${RECIPE_DIR}/../" ${SRC_DIR}/build_framework
+
+cd ${SRC_DIR}/build_framework
+cmake ${RECIPE_DIR}/../ -DCONDA_BUILD=ON \
+                        -DCMAKE_BUILD_TYPE="Release"\
+                        -DLIBRARY_LIB=$CONDA_PREFIX/lib \
+                        -DLIBRARY_INC=$CONDA_PREFIX \
+                        -DCMAKE_INSTALL_PREFIX=$PREFIX
+
+make install
+# $PYTHON setup.py install
diff --git a/recipe/conda_build_config.yaml b/recipe/conda_build_config.yaml
new file mode 100644
index 0000000..e5e168b
--- /dev/null
+++ b/recipe/conda_build_config.yaml
@@ -0,0 +1,10 @@
+python:
+  - 2.7 # [not win]
+  - 3.5
+  - 3.6
+numpy:
+  # TODO investigage, as it doesn't currently build with cvxp, requires >1.14
+  - 1.11
+  - 1.12
+  - 1.14
+  - 1.15
diff --git a/recipe/meta.yaml b/recipe/meta.yaml
new file mode 100644
index 0000000..fee8e87
--- /dev/null
+++ b/recipe/meta.yaml
@@ -0,0 +1,45 @@
+package:
+  name: ccpi-framework
+  version: {{ environ['CIL_VERSION'] }}
+
+build:
+  preserve_egg_dir: False
+  script_env:
+    - CIL_VERSION   
+  #number: 0  
+  
+test:
+  requires:
+    - python-wget
+    - cvxpy # [ unix and py36 and np115 ]
+
+  source_files:
+    - ./Wrappers/Python/test # [win]
+    - ./ccpi/Wrappers/Python/test # [not win] 
+
+  commands:
+    - python -c "import os; print ('TESTING IN THIS DIRECTORY' , os.getcwd())"
+    - python -m unittest discover Wrappers/Python/test # [win]
+    - python -m unittest discover -s ccpi/Wrappers/Python/test # [not win]
+    
+requirements:
+  build:
+    - {{ pin_compatible('numpy', max_pin='x.x') }}
+    - python
+    - setuptools
+    - cmake
+
+  run:
+    - {{ pin_compatible('numpy', max_pin='x.x') }}
+    - python
+    - numpy
+    - scipy
+    - matplotlib
+    - h5py
+    - pillow
+    - libgcc-ng # [not win]
+  
+about:
+  home: http://www.ccpi.ac.uk
+  license:  Apache 2.0 License
+  summary: 'CCPi Framework'
diff --git a/src/Core/CMakeLists.txt b/src/Core/CMakeLists.txt
new file mode 100644
index 0000000..a527289
--- /dev/null
+++ b/src/Core/CMakeLists.txt
@@ -0,0 +1,124 @@
+
+set (CMAKE_C_STANDARD 11)
+
+if(APPLE)
+    if(NOT DEFINED OPENMP_INCLUDES OR NOT DEFINED OPENMP_LIBRARIES)
+        set(OPENMP_INCLUDES "OPENMP_INCLUDES-NOTFOUND" CACHE PATH "Need to set OpenMP include dir on OSX")
+        set(OPENMP_LIBRARIES "OPENMP_LIBRARIES-NOTFOUND" CACHE PATH "Need to set OpenMP lib dir on OSX")
+    endif()
+    if(NOT EXISTS ${OPENMP_INCLUDES} OR NOT EXISTS ${OPENMP_LIBRARIES})
+        message(FATAL_ERROR "Need to set OPENMP_INCLUDES and OPENMP_LIBRARIES on OSX.")
+    endif()
+    if(CMAKE_C_COMPILER_ID MATCHES "Clang")
+        set(OpenMP_C "${CMAKE_C_COMPILER}")
+        set(OpenMP_C_FLAGS "-fopenmp=libomp -Wno-unused-command-line-argument")
+        set(OpenMP_C_LIB_NAMES "libomp" "libgomp" "libiomp5")
+        set(OpenMP_libomp_LIBRARY ${OpenMP_C_LIB_NAMES})
+        set(OpenMP_libgomp_LIBRARY ${OpenMP_C_LIB_NAMES})
+        set(OpenMP_libiomp5_LIBRARY ${OpenMP_C_LIB_NAMES})
+    endif()
+    if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+        set(OpenMP_CXX "${CMAKE_CXX_COMPILER}")
+        set(OpenMP_CXX_FLAGS "-fopenmp=libomp -Wno-unused-command-line-argument")
+        set(OpenMP_CXX_LIB_NAMES "libomp" "libgomp" "libiomp5")
+        set(OpenMP_libomp_LIBRARY ${OpenMP_CXX_LIB_NAMES})
+        set(OpenMP_libgomp_LIBRARY ${OpenMP_CXX_LIB_NAMES})
+        set(OpenMP_libiomp5_LIBRARY ${OpenMP_CXX_LIB_NAMES})
+    endif()
+endif()
+
+
+find_package(OpenMP REQUIRED)
+
+add_definitions(${OpenMP_CXX_FLAGS})
+add_definitions(${OpenMP_C_FLAGS})
+  
+  if (APPLE)
+    set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+    set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+    set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
+    set (CMAKE_SHARED_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_SHARED_LINKER_FLAGS}")
+    set (CMAKE_STATIC_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_STATIC_LINKER_FLAGS}")   
+    include_directories("${OPENMP_INCLUDES}")
+    link_directories("${OPENMP_LIBRARIES}")
+  else()
+    if(${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.9.0") 
+      set (OpenMP_EXE_LINKER_FLAGS OpenMP::OpenMP_CXX)
+      set (OpenMP_EXE_LINKER_FLAGS OpenMP::OpenMP_C)
+    else()
+      message(WARNING "Your CMake version is old. OpenMP linking flags  might be incorrect.")
+      # need to explicitly set this. Definitely for gcc, hopefully also for other systems.
+      # See https://gitlab.kitware.com/cmake/cmake/issues/15392
+      set (OpenMP_EXE_LINKER_FLAGS ${OpenMP_CXX_FLAGS})
+      set (OpenMP_EXE_LINKER_FLAGS ${OpenMP_C_FLAGS})
+    endif()
+  endif()
+
+
+#if (OPENMP_FOUND)
+#    set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+#    set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+#    set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}")
+#   set (CMAKE_SHARED_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_SHARED_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}")
+#   set (CMAKE_STATIC_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_STATIC_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}")
+#   if (UNIX)
+#    set (FLAGS "-O3 -funsigned-char -Wall  -Wl,--no-undefined -march=native")  
+#    set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}")
+#    set (CMAKE_C_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}")
+#  
+#    set (EXTRA_LIBRARIES 
+#      "gomp"
+#      "m"
+#      )
+#   endif()
+# endif()
+
+if (WIN32)
+	set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} /Ddll_EXPORTS")
+endif()
+
+
+message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}")
+message("CMAKE_C_FLAGS ${CMAKE_C_FLAGS}")
+message("CMAKE_EXE_LINKER_FLAGS ${CMAKE_EXE_LINKER_FLAGS}")
+message("CMAKE_SHARED_LINKER_FLAGS ${CMAKE_SHARED_LINKER_FLAGS}")
+message("CMAKE_STATIC_LINKER_FLAGS ${CMAKE_STATIC_LINKER_FLAGS}")
+
+
+
+add_library(cilacc SHARED ${CMAKE_CURRENT_SOURCE_DIR}/axpby.c
+                          ${CMAKE_CURRENT_SOURCE_DIR}/FiniteDifferenceLibrary.c )
+
+target_link_libraries(cilacc ${OpenMP_C_LIB_NAMES} )
+include_directories(cilacc PUBLIC 
+				  ${CMAKE_CURRENT_SOURCE_DIR}/include
+		          )
+
+## Install
+#include(GNUInstallDirs)
+#install(TARGETS cilacc
+#	RUNTIME DESTINATION bin
+#    LIBRARY DESTINATION lib
+#    ARCHIVE DESTINATION lib
+#	CONFIGURATIONS ${CMAKE_BUILD_TYPE} 
+#	)
+
+if (UNIX)
+message ("I'd install into ${CMAKE_INSTALL_PREFIX}/lib")
+install(TARGETS cilacc
+	LIBRARY DESTINATION lib
+	CONFIGURATIONS ${CMAKE_BUILD_TYPE} 
+	)
+elseif(WIN32)
+message ("I'd install into ${CMAKE_INSTALL_PREFIX} lib bin")
+  install(TARGETS cilacc 
+	RUNTIME DESTINATION bin
+	ARCHIVE DESTINATION lib
+	CONFIGURATIONS ${CMAKE_BUILD_TYPE} 
+	)
+endif()
+
+install(DIRECTORY ${PROJECT_SOURCE_DIR}/src/Core/include/
+        DESTINATION ${CMAKE_INSTALL_PREFIX}/include/cil) 
+
+
diff --git a/src/Core/FiniteDifferenceLibrary.c b/src/Core/FiniteDifferenceLibrary.c
new file mode 100644
index 0000000..f5736e2
--- /dev/null
+++ b/src/Core/FiniteDifferenceLibrary.c
@@ -0,0 +1,692 @@
+#define dll_EXPORTS = 1
+
+#include "FiniteDifferenceLibrary.h"
+
+DLL_EXPORT int openMPtest(void)
+{
+	int nThreads;
+#pragma omp parallel
+	{
+		if (omp_get_thread_num() == 0)
+		{
+			nThreads = omp_get_num_threads();
+		}
+	}
+	return nThreads;
+}
+int fdiff_direct_neumann(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, long nx, long ny, long nz, long nc)
+{
+	size_t volume = nx * ny * nz;
+
+	const float * inimage = inimagefull;
+	float * outimageX = outimageXfull;
+	float * outimageY = outimageYfull;
+	float * outimageZ = outimageZfull;
+
+	int offset1 = (nz - 1) * nx * ny; //ind to beginning of last slice
+	int offset2 = offset1 + (ny - 1) * nx; //ind to beginning of last row
+
+	long c;
+	for (c = 0; c < nc; c++)
+	{
+#pragma omp parallel
+		{
+			long ind, k, j, i;
+			float pix0;
+			//run over all and then fix boundaries
+#pragma omp for
+			for (ind = 0; ind < nx * ny * (nz - 1); ind++)
+			{
+				pix0 = -inimage[ind];
+
+				outimageX[ind] = pix0 + inimage[ind + 1];
+				outimageY[ind] = pix0 + inimage[ind + nx];
+				outimageZ[ind] = pix0 + inimage[ind + nx * ny];
+			}
+
+#pragma omp for
+			for (ind = 0; ind < nx * (ny - 1); ind++)
+			{
+				pix0 = -inimage[ind + offset1];
+
+				outimageX[ind + offset1] = pix0 + inimage[ind + offset1 + 1];
+				outimageY[ind + offset1] = pix0 + inimage[ind + offset1 + nx];
+			}
+
+#pragma omp for
+			for (ind = 0; ind < nx - 1; ind++)
+			{
+				pix0 = -inimage[ind + offset2];
+
+				outimageX[ind + offset2] = pix0 + inimage[ind + offset2 + 1];
+			}
+
+			//boundaries
+#pragma omp for
+			for (k = 0; k < nz; k++)
+			{
+				for (i = 0; i < nx; i++)
+				{
+					outimageY[(k * ny * nx) + (ny - 1) * nx + i] = 0;
+				}
+			}
+
+#pragma omp for
+			for (k = 0; k < nz; k++)
+			{
+				for (j = 0; j < ny; j++)
+				{
+					outimageX[k * ny * nx + j * nx + nx - 1] = 0;
+				}
+			}
+
+			if (nz > 1)
+			{
+#pragma omp for
+				for (ind = 0; ind < ny * nx; ind++)
+				{
+					outimageZ[nx * ny * (nz - 1) + ind] = 0;
+				}
+			}
+		}
+
+		inimage += volume;
+		outimageX += volume;
+		outimageY += volume;
+		outimageZ += volume;
+	}
+
+
+	//now the rest of the channels
+	if (nc > 1)
+	{
+		long ind;
+
+		for (c = 0; c < nc - 1; c++)
+		{
+			float * outimageC = outimageCfull + c * volume;
+			const float * inimage = inimagefull + c * volume;
+
+#pragma omp parallel for
+			for (ind = 0; ind < volume; ind++)
+			{
+				outimageC[ind] = -inimage[ind] + inimage[ind + volume];
+			}
+		}
+
+#pragma omp parallel for
+		for (ind = 0; ind < volume; ind++)
+		{
+			outimageCfull[(nc - 1) * volume + ind] = 0;
+		}
+	}
+
+	return 0;
+}
+int fdiff_direct_neumann_L21(const float *inimagefull, float *outimageL21normfull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, long nx, long ny, long nz, long nc)
+{
+	size_t volume = nx * ny * nz;
+	int z_dim = nz - 1;
+
+	const float * inimage = inimagefull;
+	float * outimageX = outimageXfull;
+	float * outimageY = outimageYfull;
+	float * outimageZ = outimageZfull;
+	float * outimageL21norm = outimageL21normfull;
+
+
+	int offset1 = (nz - 1) * nx * ny; //ind to beginning of last slice
+	int offset2 = offset1 + (ny - 1) * nx; //ind to beginning of last row
+
+	long c;
+	for (c = 0; c < nc; c++)
+	{
+#pragma omp parallel
+		{
+			long ind, k;
+			float pix0;
+			//run over all and then fix boundaries
+#pragma omp for
+			for (ind = 0; ind < nx * ny * (nz - 1); ind++)
+			{
+				pix0 = -inimage[ind];
+
+				outimageX[ind] = pix0 + inimage[ind + 1];
+				outimageY[ind] = pix0 + inimage[ind + nx];
+				outimageZ[ind] = pix0 + inimage[ind + nx * ny];
+			}
+
+#pragma omp for
+			for (ind = 0; ind < nx * (ny - 1); ind++)
+			{
+				pix0 = -inimage[ind + offset1];
+
+				outimageX[ind + offset1] = pix0 + inimage[ind + offset1 + 1];
+				outimageY[ind + offset1] = pix0 + inimage[ind + offset1 + nx];
+			}
+
+#pragma omp for
+			for (ind = 0; ind < nx - 1; ind++)
+			{
+				pix0 = -inimage[ind + offset2];
+
+				outimageX[ind + offset2] = pix0 + inimage[ind + offset2 + 1];
+			}
+
+			//boundaries
+#pragma omp for
+			for (k = 0; k < nz; k++)
+			{
+				for (int i = 0; i < nx; i++)
+				{
+					outimageY[(k * ny * nx) + (ny - 1) * nx + i] = 0;
+				}
+			}
+
+#pragma omp for
+			for (k = 0; k < nz; k++)
+			{
+				for (int j = 0; j < ny; j++)
+				{
+					outimageX[k * ny * nx + j * nx + nx - 1] = 0;
+				}
+			}
+
+			if (z_dim)
+			{
+#pragma omp for
+				for (ind = 0; ind < ny * nx; ind++)
+				{
+					outimageZ[nx * ny * (nz - 1) + ind] = 0;
+				}
+
+#pragma omp for
+				for (ind = 0; ind < volume; ind++)
+				{
+					outimageL21norm[ind] = outimageX[ind] * outimageX[ind] + outimageY[ind] * outimageY[ind] + outimageZ[ind] * outimageZ[ind];
+				}
+
+			}
+			else
+			{
+
+#pragma omp for
+				for (ind = 0; ind < volume; ind++)
+				{
+					outimageL21norm[ind] = outimageX[ind] * outimageX[ind] + outimageY[ind] * outimageY[ind];
+				}
+			}
+		}
+
+	
+		inimage += volume;
+		outimageX += volume;
+		outimageY += volume;
+		outimageZ += volume;
+		outimageL21norm += volume;
+	}
+
+
+	//now the rest of the channels
+	if (nc > 1)
+	{
+		long ind;
+
+		for (c = 0; c < nc - 1; c++)
+		{
+			float * outimageC = outimageCfull + c * volume;
+			float * outimageL21norm = outimageL21normfull + c * volume;
+			const float * inimage = inimagefull + c * volume;
+
+#pragma omp parallel for
+			for (ind = 0; ind < volume; ind++)
+			{
+				outimageC[ind] = -inimage[ind] + inimage[ind + volume];
+				outimageL21norm[ind] += outimageC[ind] * outimageC[ind];
+
+				//sqrt
+				outimageL21norm[ind] = (float)sqrt(outimageL21norm[ind]);
+			}
+		}
+
+#pragma omp parallel for
+		for (ind = 0; ind < volume; ind++)
+		{
+			outimageCfull[(nc - 1) * volume + ind] = 0;
+		}
+	}
+
+	return 0;
+}
+int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, long nx, long ny, long nz, long nc)
+{
+	size_t volume = nx * ny * nz;
+
+	const float * inimage = inimagefull;
+	float * outimageX = outimageXfull;
+	float * outimageY = outimageYfull;
+	float * outimageZ = outimageZfull;
+
+	int offset1 = (nz - 1) * nx * ny; //ind to beginning of last slice
+	int offset2 = offset1 + (ny - 1) * nx; //ind to beginning of last row
+
+
+	long c;
+	for (c = 0; c < nc; c++)
+	{
+
+#pragma omp parallel
+		{
+			long ind, k;
+			float pix0;
+			//run over all and then fix boundaries
+#pragma omp for
+			for (ind = 0; ind < nx * ny * (nz - 1); ind++)
+			{
+				pix0 = -inimage[ind];
+
+				outimageX[ind] = pix0 + inimage[ind + 1];
+				outimageY[ind] = pix0 + inimage[ind + nx];
+				outimageZ[ind] = pix0 + inimage[ind + nx * ny];
+			}
+
+#pragma omp for
+			for (ind = 0; ind < nx * (ny - 1); ind++)
+			{
+				pix0 = -inimage[ind + offset1];
+
+				outimageX[ind + offset1] = pix0 + inimage[ind + offset1 + 1];
+				outimageY[ind + offset1] = pix0 + inimage[ind + offset1 + nx];
+			}
+
+#pragma omp for
+			for (ind = 0; ind < nx - 1; ind++)
+			{
+				pix0 = -inimage[ind + offset2];
+
+				outimageX[ind + offset2] = pix0 + inimage[ind + offset2 + 1];
+			}
+
+			//boundaries
+#pragma omp for
+			for (k = 0; k < nz; k++)
+			{
+				for (int i = 0; i < nx; i++)
+				{
+					int ind1 = (k * ny * nx);
+					int ind2 = ind1 + (ny - 1) * nx;
+
+					outimageY[ind2 + i] = -inimage[ind2 + i] + inimage[ind1 + i];
+				}
+			}
+
+#pragma omp for
+			for (k = 0; k < nz; k++)
+			{
+				for (int j = 0; j < ny; j++)
+				{
+					int ind1 = k * ny * nx + j * nx;
+					int ind2 = ind1 + nx - 1;
+
+					outimageX[ind2] = -inimage[ind2] + inimage[ind1];
+				}
+			}
+
+
+			if (nz > 1)
+			{
+#pragma omp for
+				for (ind = 0; ind < ny * nx; ind++)
+				{
+					outimageZ[nx * ny * (nz - 1) + ind] = -inimage[nx * ny * (nz - 1) + ind] + inimage[ind];
+				}
+			}
+		}
+
+		inimage += volume;
+		outimageX += volume;
+		outimageY += volume;
+		outimageZ += volume;
+	}
+
+	//now the rest of the channels
+	if (nc > 1)
+	{
+		long ind;
+
+		for (c = 0; c < nc - 1; c++)
+		{
+			float * outimageC = outimageCfull + c * volume;
+			const float * inimage = inimagefull + c * volume;
+
+#pragma omp parallel for
+			for (ind = 0; ind < volume; ind++)
+			{
+				outimageC[ind] = -inimage[ind] + inimage[ind + volume];
+			}
+		}
+
+#pragma omp parallel for
+		for (ind = 0; ind < volume; ind++)
+		{
+			outimageCfull[(nc - 1) * volume + ind] = -inimagefull[(nc - 1) * volume + ind] + inimagefull[ind];
+		}
+	}
+
+	return 0;
+}
+int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const float *inimageYfull, const float *inimageZfull, const float *inimageCfull, long nx, long ny, long nz, long nc)
+{
+	//runs over full data in x, y, z. then corrects elements for bounday conditions and sums
+	size_t volume = nx * ny * nz;
+
+	//assumes nx and ny > 1
+	int z_dim = nz - 1;
+
+	float * outimage = outimagefull;
+	const float * inimageX = inimageXfull;
+	const float * inimageY = inimageYfull;
+	const float * inimageZ = inimageZfull;
+
+	float * tempX = (float*)malloc(volume * sizeof(float));
+	float * tempY = (float*)malloc(volume * sizeof(float));
+	float * tempZ;
+	
+	if(z_dim)
+	{
+		tempZ = (float*)malloc(volume * sizeof(float));
+	}
+
+	long c;
+	for (c = 0; c < nc; c++) //just calculating x, y and z in each channel here
+	{
+#pragma omp parallel
+		{
+			long ind, k;
+
+#pragma omp for
+			for (ind = 1; ind < nx * ny * nz; ind++)
+			{
+				tempX[ind] = -inimageX[ind] + inimageX[ind - 1];
+			}
+#pragma omp for
+			for (ind = nx; ind < nx * ny * nz; ind++)
+			{
+				tempY[ind] = -inimageY[ind] + inimageY[ind - nx];
+
+			}
+
+			//boundaries
+#pragma omp for
+			for (k = 0; k < nz; k++)
+			{
+				for (int j = 0; j < ny; j++)
+				{
+					tempX[k * ny * nx + j * nx] = -inimageX[k * ny * nx + j * nx];
+					tempX[k * ny * nx + j * nx + nx - 1] = inimageX[k * ny * nx + j * nx + nx - 2];
+				}
+			}
+#pragma omp for
+			for (k = 0; k < nz; k++)
+			{
+				for (int i = 0; i < nx; i++)
+				{
+					tempY[(k * ny * nx) + i] = -inimageY[(k * ny * nx) + i];
+					tempY[(k * ny * nx) + nx * (ny - 1) + i] = inimageY[(k * ny * nx) + nx * (ny - 2) + i];
+				}
+			}
+
+			if (z_dim)
+			{
+#pragma omp for
+				for (ind = nx * ny; ind < nx * ny * nz; ind++)
+				{
+					tempZ[ind] = -inimageZ[ind] + inimageZ[ind - nx * ny];
+				}				
+#pragma omp for
+				for (ind = 0; ind < ny * nx; ind++)
+				{
+					tempZ[ind] = -inimageZ[ind];
+					tempZ[nx * ny * (nz - 1) + ind] = inimageZ[nx * ny * (nz - 2) + ind];
+				}
+#pragma omp for
+				for (ind = 0; ind < volume; ind++)
+				{
+					outimage[ind] = tempX[ind] + tempY[ind] + tempZ[ind];
+				}
+
+			}
+			else
+			{
+#pragma omp for
+				for (ind = 0; ind < volume; ind++)
+				{
+					outimage[ind] = tempX[ind] + tempY[ind];
+				}
+			}
+
+		}
+
+		outimage += volume;
+		inimageX += volume;
+		inimageY += volume;
+		inimageZ += volume;
+	}
+	free(tempX);
+	free(tempY);
+
+	if(z_dim)
+		free(tempZ);
+
+
+	//	//now the rest of the channels
+	if (nc > 1)
+	{
+		long ind;
+
+		for (c = 1; c < nc - 1; c++)
+		{
+
+#pragma omp parallel for
+			for (ind = 0; ind < volume; ind++)
+			{
+				outimagefull[ind + c * volume] += -inimageCfull[ind + c * volume] + inimageCfull[ind + (c - 1) * volume];
+			}
+		}
+
+#pragma omp parallel for
+		for (ind = 0; ind < volume; ind++)
+		{
+			outimagefull[ind] += -inimageCfull[ind];
+			outimagefull[(nc - 1) * volume + ind] += inimageCfull[(nc - 2) * volume + ind];
+		}
+
+	}
+
+	return 0;
+}
+int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const float *inimageYfull, const float *inimageZfull, const float *inimageCfull, long nx, long ny, long nz, long nc)
+{
+	//runs over full data in x, y, z. then correctects elements for bounday conditions and sums
+	size_t volume = nx * ny * nz;
+
+	//assumes nx and ny > 1
+	int z_dim = nz - 1;
+
+	float * outimage = outimagefull;
+	const float * inimageX = inimageXfull;
+	const float * inimageY = inimageYfull;
+	const float * inimageZ = inimageZfull;
+
+	float * tempX = (float*)malloc(volume * sizeof(float));
+	float * tempY = (float*)malloc(volume * sizeof(float));
+	float * tempZ;
+
+	if (z_dim)
+	{
+		tempZ = (float*)malloc(volume * sizeof(float));
+	}
+
+	long c;
+	for (c = 0; c < nc; c++) //just calculating x, y and z in each channel here
+	{
+#pragma omp parallel
+		{
+			long ind, k;
+
+			//run over all and then fix boundaries
+#pragma omp for
+			for (ind = 1; ind < volume; ind++)
+			{
+				tempX[ind] = -inimageX[ind] + inimageX[ind - 1];
+			}
+#pragma omp for
+			for (ind = nx; ind < volume; ind++)
+			{
+				tempY[ind] = -inimageY[ind] + inimageY[ind - nx];
+			}
+
+			//boundaries
+#pragma omp for
+			for (k = 0; k < nz; k++)
+			{
+				for (int i = 0; i < nx; i++)
+				{
+					tempY[(k * ny * nx) + i] = -inimageY[(k * ny * nx) + i] + inimageY[(k * ny * nx) + nx * (ny - 1) + i];
+				}
+			}
+#pragma omp for
+			for (k = 0; k < nz; k++)
+			{
+				for (int j = 0; j < ny; j++)
+				{
+					tempX[k * ny * nx + j * nx] = -inimageX[k * ny * nx + j * nx] + inimageX[k * ny * nx + j * nx + nx - 1];
+				}
+			}
+
+			if (z_dim)
+			{
+
+#pragma omp for
+				for (ind = nx * ny; ind < nx * ny * nz; ind++)
+				{
+					tempZ[ind] = -inimageZ[ind] + inimageZ[ind - nx * ny];
+				}
+#pragma omp for
+				for (ind = 0; ind < ny * nx; ind++)
+				{
+					tempZ[ind] = -inimageZ[ind] + inimageZ[nx * ny * (nz - 1) + ind];
+				}
+
+#pragma omp for
+				for (ind = 0; ind < volume; ind++)
+				{
+					outimage[ind] = tempX[ind] + tempY[ind] + tempZ[ind];
+				}
+
+			}
+			else
+			{
+#pragma omp for
+				for (ind = 0; ind < volume; ind++)
+				{
+					outimage[ind] = tempX[ind] + tempY[ind];
+				}
+			}
+
+		}
+
+		outimage += volume;
+		inimageX += volume;
+		inimageY += volume;
+		inimageZ += volume;
+	}
+	free(tempX);
+	free(tempY);
+
+	if(z_dim)
+		free(tempZ);
+
+	//now the rest of the channels
+	if (nc > 1)
+	{
+		long ind;
+
+		for (c = 1; c < nc; c++)
+		{
+
+#pragma omp parallel for
+			for (ind = 0; ind < volume; ind++)
+			{
+				outimagefull[ind + c * volume] += -inimageCfull[ind + c * volume] + inimageCfull[ind + (c - 1) * volume];
+			}
+		}
+
+#pragma omp parallel for
+		for (ind = 0; ind < volume; ind++)
+		{
+			outimagefull[ind] += -inimageCfull[ind] + inimageCfull[(nc - 1) * volume + ind];
+		}
+	}
+
+	return 0;
+}
+
+
+DLL_EXPORT int fdiff4D(float *imagefull, float *gradCfull, float *gradZfull, float *gradYfull, float *gradXfull, long nc, long nz, long ny, long nx, int boundary, int direction)
+{
+	if (boundary)
+	{
+		if (direction)
+			fdiff_direct_periodic(imagefull, gradXfull, gradYfull, gradZfull, gradCfull, nx, ny, nz, nc);
+		else
+			fdiff_adjoint_periodic(imagefull, gradXfull, gradYfull, gradZfull, gradCfull, nx, ny, nz, nc);
+	}
+	else
+		{
+		if (direction)
+			fdiff_direct_neumann(imagefull, gradXfull, gradYfull, gradZfull, gradCfull, nx, ny, nz, nc);
+		else
+			fdiff_adjoint_neumann(imagefull, gradXfull, gradYfull, gradZfull, gradCfull, nx, ny, nz, nc);
+	}
+	
+	return 0;
+}
+DLL_EXPORT int fdiff3D(float *imagefull, float *gradZfull, float *gradYfull, float *gradXfull, long nz, long ny, long nx, int boundary, int direction)
+{
+	if (boundary)
+	{
+		if (direction)
+			fdiff_direct_periodic(imagefull, gradXfull, gradYfull, gradZfull, NULL, nx, ny, nz, 1);
+		else
+			fdiff_adjoint_periodic(imagefull, gradXfull, gradYfull, gradZfull, NULL, nx, ny, nz, 1);
+	}
+	else
+	{
+		if (direction)
+			fdiff_direct_neumann(imagefull, gradXfull, gradYfull, gradZfull, NULL, nx, ny, nz, 1);
+		else
+			fdiff_adjoint_neumann(imagefull, gradXfull, gradYfull, gradZfull, NULL, nx, ny, nz, 1);
+	}
+
+	return 0;
+}
+DLL_EXPORT int fdiff2D(float *imagefull, float *gradYfull, float *gradXfull, long ny, long nx, int boundary, int direction)
+{
+	if (boundary)
+	{
+		if (direction)
+			fdiff_direct_periodic(imagefull, gradXfull, gradYfull, NULL, NULL, nx, ny, 1, 1);
+		else
+			fdiff_adjoint_periodic(imagefull, gradXfull, gradYfull, NULL, NULL, nx, ny, 1, 1);
+	}
+	else
+	{
+		if (direction)
+			fdiff_direct_neumann(imagefull, gradXfull, gradYfull, NULL, NULL, nx, ny, 1, 1);
+		else
+			fdiff_adjoint_neumann(imagefull, gradXfull, gradYfull, NULL, NULL, nx, ny, 1, 1);
+	}
+
+	return 0;
+}
diff --git a/src/Core/axpby.c b/src/Core/axpby.c
new file mode 100755
index 0000000..c4d162d
--- /dev/null
+++ b/src/Core/axpby.c
@@ -0,0 +1,109 @@
+#include "axpby.h"
+
+
+DLL_EXPORT int padd(float * x, float * y, float * out, long size){
+    long i = 0;
+#pragma omp parallel for
+    for (i=0; i < size; i++)
+    {
+        *(out + i ) = *(x + i) + *(y+i); 
+    }
+    return 0;
+}
+
+DLL_EXPORT int psubtract(float * x, float * y, float * out, long size){
+    long i = 0;
+#pragma omp parallel
+{
+//#pragma omp single
+//{
+//		printf("current number of threads %d\n", omp_get_num_threads());
+//}
+#pragma omp for
+    for (i=0; i < size; i++)
+    {
+        *(out + i ) = *(x + i) - *(y+i); 
+    }
+}
+    return 0;
+    
+}
+
+DLL_EXPORT int pmultiply(float * x, float * y, float * out, long size){
+    long i = 0;
+#pragma omp parallel for
+    for (i=0; i < size; i++)
+    {
+        *(out + i ) = *(x + i) * *(y+i); 
+    }
+    return 0;
+}
+
+DLL_EXPORT int pdivide(float * x, float * y, float * out, long size, float default_value)
+{
+    long i = 0;
+#pragma omp parallel for
+    for (i=0; i < size; i++)
+    {
+        *(out + i ) = *(y+i) ? *(x + i) / *(y+i) : default_value;
+    }
+    return 0;
+}
+DLL_EXPORT int ppower(float * x, float * y, float * out, long size){
+    long i = 0;
+#pragma omp parallel for
+    for (i=0; i < size; i++)
+    {
+        *(out + i ) = (float)pow(*(x + i) , *(y+i)) ; 
+    }
+    return 0;
+}
+
+DLL_EXPORT int pminimum(float * x, float * y, float * out, long size){
+    long i = 0;
+#pragma omp parallel for
+    for (i=0; i < size; i++)
+    {
+        *(out + i ) = *(y+i) > (*x+i) ? *(x + i) :  *(y+i); 
+    }
+    return 0;
+}
+
+DLL_EXPORT int pmaximum(float * x, float * y, float * out, long size) {
+	long i = 0;
+#pragma omp parallel for
+	for (i = 0; i < size; i++)
+	{
+		*(out + i) = *(y + i) < (*x + i) ? *(x + i) : *(y + i);
+	}
+	return 0;
+}
+
+
+DLL_EXPORT int saxpby(float * x, float * y, float * out, float a, float b, long size){
+    long i = 0;
+#pragma omp parallel
+{
+#pragma omp for
+    for (i=0; i < size; i++)
+    {
+        *(out + i ) = a * ( *(x + i) ) + b * ( *(y + i) ); 
+    }
+}
+    return 0;
+    
+}
+
+DLL_EXPORT int daxpby(double * x, double * y, double * out, double a, double b, long size) {
+	long i = 0;
+#pragma omp parallel
+	{
+#pragma omp for
+		for (i = 0; i < size; i++)
+		{
+			*(out + i) = a * (*(x + i)) + b * (*(y + i));
+		}
+	}
+	return 0;
+
+}
\ No newline at end of file
diff --git a/src/Core/include/FiniteDifferenceLibrary.h b/src/Core/include/FiniteDifferenceLibrary.h
new file mode 100644
index 0000000..6e426af
--- /dev/null
+++ b/src/Core/include/FiniteDifferenceLibrary.h
@@ -0,0 +1,6 @@
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "omp.h"
+//#include "ipp.h"
+#include "dll_export.h"
\ No newline at end of file
diff --git a/src/Core/include/axpby.h b/src/Core/include/axpby.h
new file mode 100644
index 0000000..2849547
--- /dev/null
+++ b/src/Core/include/axpby.h
@@ -0,0 +1,17 @@
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "omp.h"
+#include "dll_export.h"
+
+
+DLL_EXPORT int padd(float * x, float * y, float * out, long size);
+DLL_EXPORT int psubtract(float * x, float * y, float * out, long size);
+DLL_EXPORT int pmultiply(float * x, float * y, float * out, long size);
+DLL_EXPORT int pdivide(float * x, float * y, float * out, long size, float default_value);
+DLL_EXPORT int ppower(float * x, float * y, float * out, long size);
+DLL_EXPORT int pminimum(float * x, float * y, float * out, long size);
+DLL_EXPORT int pmaximum(float * x, float * y, float * out, long size);
+
+DLL_EXPORT int saxpby(float * x, float * y, float * out, float a, float b, long size);
+DLL_EXPORT int daxpby(double * x, double * y, double * out, double a, double b, long size);
diff --git a/src/Core/include/dll_export.h b/src/Core/include/dll_export.h
new file mode 100755
index 0000000..6b816c3
--- /dev/null
+++ b/src/Core/include/dll_export.h
@@ -0,0 +1,20 @@
+#pragma once
+#ifndef DLLEXPORT_H
+#define DLLEXPORT_H
+
+#if defined(_WIN32) || defined(__WIN32__)
+#if defined(dll_EXPORTS)  // add by CMake 
+#define  DLL_EXPORT __declspec(dllexport)
+#define EXPIMP_TEMPLATE
+#else
+#define  DLL_EXPORT __declspec(dllimport)
+#define EXPIMP_TEMPLATE extern
+#endif 
+#elif defined(linux) || defined(__linux) || defined(__APPLE__)
+#define DLL_EXPORT
+#ifndef __cdecl
+#define __cdecl
+#endif
+#endif
+
+#endif