Merge pull request #14 from joeljonsson/develop

update master

Author: cheesema

.gitmodules

@@ -1,7 +1,8 @@
 [submodule "LibAPR"]
 	path = LibAPR
 	url = https://github.com/AdaptiveParticles/LibAPR.git
-	branch = develop
+	branch = develop_joel
 [submodule "pybind11"]
 	path = pybind11
 	url = https://github.com/pybind/pybind11.git
+	branch = master

CMakeLists.txt

@@ -6,6 +6,7 @@ set(CMAKE_CXX_STANDARD 14)
 
 # APR build options:
 option(PYAPR_USE_OPENMP "Use OpenMP for PyAPR functionality" ON)
+option(PYAPR_USE_CUDA "Use CUDA for PyAPR functionality" OFF)
 
 ##------------- Required libraries -------------##
 find_package(HDF5 REQUIRED)
@@ -23,16 +24,16 @@ set(APR_TESTS OFF CACHE BOOL "" FORCE)
 set(APR_BUILD_STATIC_LIB ON CACHE BOOL "" FORCE)
 set(APR_BUILD_SHARED_LIB OFF CACHE BOOL "" FORCE)
 set(APR_BUILD_EXAMPLES OFF CACHE BOOL "" FORCE)
-set(APR_USE_CUDA OFF CACHE BOOL "" FORCE)
 set(APR_BUILD_PYTHON_WRAPPERS OFF CACHE BOOL "" FORCE)
 set(APR_USE_LIBTIFF ON CACHE BOOL "" FORCE)
-set(APR_USE_OPENMP ON CACHE BOOL "" FORCE)
+set(APR_USE_OPENMP ${PYAPR_USE_OPENMP} CACHE BOOL "" FORCE)
+set(APR_USE_CUDA ${PYAPR_USE_CUDA} CACHE BOOL "" FORCE)
 
 include_directories(${CMAKE_CURRENT_BINARY_DIR} ${CMAKE_CURRENT_BINARY_DIR}/LibAPR ${HDF5_INCLUDE_DIRS} ${ZLIB_INCLUDE_DIRS} ${TIFF_INCLUDE_DIR})
 
 add_subdirectory("LibAPR")
 include_directories(LibAPR/src)
-
+include_directories("LibAPR/external/glm")
 
 ###############################################################################
 # Handle OpenMP
@@ -41,11 +42,26 @@ find_package(OpenMP)
 if(NOT OPENMP_FOUND OR NOT PYAPR_USE_OPENMP)
     message(WARNING "OpenMP support not found or disabled with current compiler. While PyAPR can compile like this, performance might not be optimal. Please see README.md for instructions.")
 else()
-    set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DPYAPR_HAVE_OPENMP ${OpenMP_C_FLAGS}")
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DPYAPR_HAVE_OPENMP ${OpenMP_CXX_FLAGS}")
+    set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DPYAPR_HAVE_OPENMP -DHAVE_OPENMP ${OpenMP_C_FLAGS}")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DPYAPR_HAVE_OPENMP -DHAVE_OPENMP ${OpenMP_CXX_FLAGS}")
 endif()
 
 
+###############################################################################
+# Handle CUDA
+###############################################################################
+if(PYAPR_USE_CUDA)
+    message(STATUS "APR: Building CUDA for PYAPR")
+    set(CMAKE_CUDA_STANDARD 14)
+    set(CMAKE_CUDA_RUNTIME_LIBRARY "Static")
+    set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} --default-stream per-thread -lineinfo -Xptxas -O3,-v --use_fast_math -DPYAPR_USE_CUDA -DAPR_USE_CUDA")
+
+    enable_language(CUDA)
+    set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DPYAPR_USE_CUDA -DAPR_USE_CUDA")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DPYAPR_USE_CUDA -DAPR_USE_CUDA")
+    include_directories(${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
+
+endif()
 
 ###############################################################################
 # Configure compiler options

LibAPR

@@ -0,0 +1,74 @@
+import pyapr
+import numpy as np
+from skimage import io as skio
+
+
+def main():
+
+    # Read in an image
+    io_int = pyapr.filegui.InteractiveIO()
+    fpath = io_int.get_tiff_file_name()
+    img = skio.imread(fpath).astype(np.uint16)
+
+    # Initialize objects
+    apr = pyapr.APR()
+    parts = pyapr.ShortParticles()
+    par = pyapr.APRParameters()
+    converter = pyapr.converter.ShortConverter()
+
+    # Set some parameters
+    par.auto_parameters = False
+    par.rel_error = 0.1
+    par.Ip_th = 10
+    par.grad_th = 50
+    par.gradient_smoothing = 2
+    par.sigma_th = 100
+    par.sigma_th_max = 50
+    converter.set_parameters(par)
+    converter.set_verbose(False)
+
+    # Compute APR and sample particle values
+    converter.get_apr(apr, img)
+
+    # Compute and display the computational ratio
+    numParts = apr.total_number_particles()
+    numPix = img.size
+    CR = numPix / numParts
+
+    print('input image size: {} pixels, APR size: {} particles --> Computational Ratio: {}'.format(numPix, numParts, CR))
+
+    # Sample particle intensities
+    parts.sample_image(apr, img)
+
+    # Save the APR to file
+    fpath_apr = io_int.save_apr_file_name()
+
+    # Initialize APRFile for I/O
+    aprfile = pyapr.io.APRFile()
+    aprfile.set_read_write_tree(True)
+
+    # Write APR and particles to file
+    aprfile.open(fpath_apr, 'WRITE')
+    aprfile.write_apr(apr)
+    aprfile.write_particles('particles', parts)
+    aprfile.close()
+
+    # Read the newly written file
+
+    # Initialize objects for reading in data
+    apr2 = pyapr.APR()
+    parts2 = pyapr.ShortParticles()
+
+    # Read from APR file
+    aprfile.open(fpath_apr, 'READ')
+    aprfile.read_apr(apr2)
+    aprfile.read_particles(apr2, 'particles', parts2)
+    aprfile.close()
+
+    # Reconstruct pixel image
+    tmp = pyapr.numerics.reconstruction.recon_pc(apr, parts)
+    recon = np.array(tmp, copy=False)
+
+
+if __name__ == '__main__':
+    main()

demo/apr_io_demo.py

@@ -0,0 +1,79 @@
+import pyapr
+import numpy as np
+from time import time
+
+
+def main():
+
+    io_int = pyapr.filegui.InteractiveIO()
+    fpath_apr = io_int.get_apr_file_name()  # get APR file path from gui
+
+    aprfile = pyapr.io.APRFile()
+    aprfile.set_read_write_tree(True)
+
+    # Initialize APR and particle objects
+    parts = pyapr.ShortParticles()
+    apr = pyapr.APR()
+
+    # Read from APR file
+    aprfile.open(fpath_apr, 'READ')
+    aprfile.read_apr(apr)
+    aprfile.read_particles(apr, 'particles', parts)
+    aprfile.close()
+
+    start = time()
+
+    org_dims = apr.org_dims()  # (Ny, Nx, Nz)
+    py_recon = np.empty((org_dims[2], org_dims[1], org_dims[0]), dtype=np.uint16)
+    max_level = apr.level_max()
+
+    apr_it = apr.iterator()
+
+    # loop over levels up to level_max-1
+    for level in range(apr_it.level_min(), apr_it.level_max()):
+
+        step_size = 2 ** (max_level - level)
+
+        for z in range(apr_it.z_num(level)):
+            for x in range(apr_it.x_num(level)):
+                for idx in range(apr_it.begin(level, z, x), apr_it.end()):
+                    y = apr_it.y(idx)  # this is slow
+
+                    y_start = y * step_size
+                    x_start = x * step_size
+                    z_start = z * step_size
+
+                    y_end = min(y_start+step_size, py_recon.shape[2])
+                    x_end = min(x_start+step_size, py_recon.shape[1])
+                    z_end = min(z_start+step_size, py_recon.shape[0])
+
+                    py_recon[z_start:z_end, x_start:x_end, y_start:y_end] = parts[idx]
+
+    # particles at the maximum level coincide with pixels
+    level = max_level
+    for z in range(apr_it.z_num(level)):
+        for x in range(apr_it.x_num(level)):
+            for idx in range(apr_it.begin(level, z, x), apr_it.end()):
+                py_recon[z, x, apr_it.y(idx)] = parts[idx]
+
+    py_time = time()-start
+    print('python reconstruction took {} seconds'.format(py_time))
+
+    start = time()
+    tmp = pyapr.numerics.reconstruction.recon_pc(apr, parts)
+    cpp_recon = np.array(tmp, copy=False)
+    cpp_time = time()-start
+    print('c++ reconstruction took {} seconds'.format(cpp_time))
+    print('c++ was {} times faster'.format(py_time / cpp_time))
+
+    zm = min(org_dims[2], 128)
+    xm = min(org_dims[1], 128)
+    ym = min(org_dims[0], 128)
+
+    success = np.allclose(py_recon[:zm, :xm, :ym], cpp_recon[:zm, :xm, :ym])
+    if not success:
+        print('Python and C++ reconstructions give different results...')
+
+
+if __name__ == '__main__':
+    main()

demo/apr_iteration_demo.py

@@ -0,0 +1,50 @@
+import pyapr
+
+
+def main():
+
+    io_int = pyapr.filegui.InteractiveIO()
+    fpath_apr = io_int.get_apr_file_name()  # get APR file path from gui
+
+    # Initialize objects
+    apr = pyapr.APR()
+    parts = pyapr.ShortParticles()
+
+    # Initialize APRFile for I/O
+    aprfile = pyapr.io.APRFile()
+    aprfile.set_read_write_tree(True)
+
+    # Read APR and particles from file
+    aprfile.open(fpath_apr, 'READ')
+    aprfile.read_apr(apr)
+    aprfile.read_particles(apr, 'particles', parts)
+
+    original_file_size = aprfile.current_file_size_MB()
+    aprfile.close()
+
+    # Interactive WNL compression
+    pyapr.viewer.interactive_compression(apr, parts)
+
+    # Write compressed APR to file
+    fpath_apr_save = io_int.save_apr_file_name()  # get file path from gui
+
+    aprfile.open(fpath_apr_save, 'WRITE')
+    aprfile.write_apr(apr)  # need to write the unchanged apr structure
+    aprfile.write_particles('particles', parts)
+
+    compressed_file_size = aprfile.current_file_size_MB()
+    aprfile.close()
+
+    # Uncompressed pixel image size in MB
+    original_image_size = 2e-6 * apr.x_num(apr.level_max()) * apr.y_num(apr.level_max()) * apr.z_num(apr.level_max())
+
+    print("Original File Size: {:7.2f} MB".format(original_file_size))
+    print("Lossy Compressed File Size: {:7.2f} MB".format(compressed_file_size))
+
+    # compare uncompressed pixel image size to compressed APR file sizes
+    print("Original Memory Compression Ratio: {:7.2f} ".format(original_image_size/original_file_size))
+    print("Lossy Memory Compression Ratio: {:7.2f} ".format(original_image_size/compressed_file_size))
+
+
+if __name__ == '__main__':
+    main()