merged develop into developPython

Author: joeljonsson

CMakeLists.txt

@@ -12,6 +12,7 @@ option(APR_INSTALL "Install APR library" OFF)
 option(APR_BUILD_SHARED_LIB "Builds shared library" ON)
 option(APR_BUILD_STATIC_LIB "Builds static library" OFF)
 option(APR_BUILD_EXAMPLES "Build APR examples" OFF)
+option(APR_USE_LIBTIFF "Use LibTIFF" ON)
 option(APR_TESTS "Build APR tests" OFF)
 option(APR_PREFER_EXTERNAL_GTEST "When found, use the installed GTEST libs instead of included sources" OFF)
 option(APR_PREFER_EXTERNAL_BLOSC "When found, use the installed BLOSC libs instead of included sources" OFF)
@@ -53,7 +54,13 @@ message("Configuring for APR version: " ${APR_VERSION_STRING})
 ###############################################################################
 find_package(HDF5 REQUIRED)
 find_package(ZLIB REQUIRED)
-find_package(TIFF REQUIRED)
+if(APR_USE_LIBTIFF)
+    find_package(TIFF)
+    set(CMAKE_C_FLAGS "${CMAKE_CFLAGS} -DHAVE_LIBTIFF")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXXFLAGS} -DHAVE_LIBTIFF")
+else()
+    message(WARNING "LibTIFF support disable, this disables TIFF writing functionality.")
+endif()
 
 # Handle OpenMP
 find_package(OpenMP)
@@ -184,7 +191,11 @@ if(APR_BUILD_SHARED_LIB)
         if (APPLE)
             target_link_libraries(${SHARED_TARGET_NAME} PRIVATE -Wl,-force_load,$<TARGET_FILE:blosc_static>)
         else()
-            target_link_libraries(${SHARED_TARGET_NAME} PRIVATE -Wl,--whole-archive $<TARGET_FILE:blosc_static> -Wl,--no-whole-archive)
+            if(WIN32)
+                target_link_libraries(${SHARED_TARGET_NAME} PRIVATE blosc_static)
+            else()
+                target_link_libraries(${SHARED_TARGET_NAME} PRIVATE -Wl,--whole-archive $<TARGET_FILE:blosc_static> -Wl,--no-whole-archive)
+            endif()
         endif()
     endif()
 endif()

cmake/AddStaticLibs.cmake

@@ -1,6 +1,8 @@
 # Adds staticLib to provided (existing) static target. Useful for merging all dependencies to
 # one fat static lib.
 # Use: addStaticLibs(libStatic someStaticLibToMerge [evenMoreStaticLibsIfNeeded])
+cmake_policy(SET CMP0026 OLD)
+
 function(addStaticLibs outLibTarget)
     get_target_property(libtype ${outLibTarget} TYPE)
     if(NOT libtype STREQUAL "STATIC_LIBRARY")
@@ -17,7 +19,8 @@ function(addStaticLibs outLibTarget)
 		if(NOT libtype STREQUAL "STATIC_LIBRARY")
 			message(FATAL_ERROR "[${lib}] is not a static lib!")
 		endif()
-		list(APPEND filesToMerge $<TARGET_FILE:${lib}>)
+		get_target_property(myLib ${lib} LOCATION)
+		set(filesToMerge ${filesTomerge} ${myLib})
 	endforeach()
 
     set(outLibFile $<TARGET_FILE:${outLibTarget}>)

examples/Example_get_apr.cpp

@@ -26,6 +26,11 @@ Advanced (Direct) Settings:
 -rel_error rel_error_value (Reasonable ranges are from .08-.15), Default: 0.1
 -normalize_input (flag that will rescale the input from the input data range to 80% of the output data type range, useful for float scaled datasets)
 -compress_level (the IO uses BLOSC for lossless compression of the APR, this can be set from 1-9, where higher increases the compression level. Note, this can come at a significant time increase.)
+-compress_type (Default: 0, loss-less compression of partilce intensities, (1,2) WNL (Balázs et al. 2017) - approach compression applied to particles (1 = without prediction, 2 = with)
+
+-neighborhood_optimization_off turns off the neighborhood opetimization (This results in boundary Particle Cells also being increased in resolution after the Pulling Scheme step)
+-output_steps Writes tiff images of the individual steps (gradient magnitude, local intensity scale, and final level of the APR calculation).
+
 )";
 
 #include <algorithm>
@@ -43,31 +48,31 @@ int main(int argc, char **argv) {
     //the apr datastructure
     APR<uint16_t> apr;
 
-    APRConverter<uint16_t> apr_converter;
-
     //read in the command line options into the parameters file
-    apr_converter.par.Ip_th = options.Ip_th;
-    apr_converter.par.rel_error = options.rel_error;
-    apr_converter.par.lambda = options.lambda;
-    apr_converter.par.mask_file = options.mask_file;
-    apr_converter.par.min_signal = options.min_signal;
-    apr_converter.par.SNR_min = options.SNR_min;
-    apr_converter.par.normalized_input = options.normalize_input;
+    apr.parameters.Ip_th = options.Ip_th;
+    apr.parameters.rel_error = options.rel_error;
+    apr.parameters.lambda = options.lambda;
+    apr.parameters.mask_file = options.mask_file;
+    apr.parameters.min_signal = options.min_signal;
+    apr.parameters.SNR_min = options.SNR_min;
+    apr.parameters.normalized_input = options.normalize_input;
+    apr.parameters.neighborhood_optimization = options.neighborhood_optimization;
+    apr.parameters.output_steps = options.output_steps;
 
     //where things are
-    apr_converter.par.input_image_name = options.input;
-    apr_converter.par.input_dir = options.directory;
-    apr_converter.par.name = options.output;
-    apr_converter.par.output_dir = options.output_dir;
+    apr.parameters.input_image_name = options.input;
+    apr.parameters.input_dir = options.directory;
+    apr.parameters.name = options.output;
+    apr.parameters.output_dir = options.output_dir;
 
-    apr_converter.fine_grained_timer.verbose_flag = false;
-    apr_converter.method_timer.verbose_flag = false;
-    apr_converter.computation_timer.verbose_flag = false;
-    apr_converter.allocation_timer.verbose_flag = false;
-    apr_converter.total_timer.verbose_flag = true;
+    apr.apr_converter.fine_grained_timer.verbose_flag = false;
+    apr.apr_converter.method_timer.verbose_flag = false;
+    apr.apr_converter.computation_timer.verbose_flag = false;
+    apr.apr_converter.allocation_timer.verbose_flag = false;
+    apr.apr_converter.total_timer.verbose_flag = true;
 
     //Gets the APR
-    if(apr_converter.get_apr(apr)){
+    if(apr.get_apr()){
 
         //Below is IO and outputting of the Implied Resolution Function through the Particle Cell level.
 
@@ -79,18 +84,6 @@ int main(int argc, char **argv) {
 
         timer.verbose_flag = true;
 
-        PixelData<uint16_t> level;
-
-        apr.interp_depth_ds(level);
-
-        std::cout << std::endl;
-
-        std::cout << "Saving down-sampled Particle Cell level as tiff image" << std::endl;
-
-        std::string output_path = save_loc + file_name + "_level.tif";
-        //write output as tiff
-        TiffUtils::saveMeshAsTiff(output_path, level);
-
         std::cout << std::endl;
         float original_pixel_image_size = (2.0f*apr.orginal_dimensions(0)*apr.orginal_dimensions(1)*apr.orginal_dimensions(2))/(1000000.0);
         std::cout << "Original image size: " << original_pixel_image_size << " MB" << std::endl;
@@ -104,6 +97,8 @@ int main(int argc, char **argv) {
         unsigned int blosc_comp_level = options.compress_level;
         unsigned int blosc_shuffle = 1;
 
+        apr.apr_compress.set_compression_type(options.compress_type);
+
         //write the APR to hdf5 file
         FileSizeInfo fileSizeInfo = apr.write_apr(save_loc,file_name,blosc_comp_type,blosc_comp_level,blosc_shuffle);
         float apr_file_size = fileSizeInfo.total_file_size;
@@ -117,6 +112,20 @@ int main(int argc, char **argv) {
         std::cout << "Lossy Compression Ratio: " << original_pixel_image_size/apr_file_size << std::endl;
         std::cout << std::endl;
 
+        if(options.output_steps) {
+
+            PixelData<uint16_t> level;
+
+            apr.interp_depth_ds(level);
+
+            std::cout << std::endl;
+
+            std::cout << "Saving down-sampled Particle Cell level as tiff image" << std::endl;
+
+            std::string output_path = save_loc + file_name + "_level.tif";
+            //write output as tiff
+            TiffUtils::saveMeshAsTiff(output_path, level);
+        }
 
         } else {
         std::cout << "Oops, something went wrong. APR not computed :(." << std::endl;
@@ -220,14 +229,25 @@ cmdLineOptions read_command_line_options(int argc, char **argv){
         result.compress_level = (unsigned int)std::stoi(std::string(get_command_option(argv, argv + argc, "-compress_level")));
     }
 
+    if(command_option_exists(argv, argv + argc, "-compress_type"))
+    {
+        result.compress_type = (unsigned int)std::stoi(std::string(get_command_option(argv, argv + argc, "-compress_type")));
+    }
+
     if(command_option_exists(argv, argv + argc, "-normalize_input"))
     {
         result.normalize_input = true;
     }
 
-    if(command_option_exists(argv, argv + argc, "-store_delta"))
+    if(command_option_exists(argv, argv + argc, "-neighborhood_optimization_off"))
+    {
+        result.neighborhood_optimization = false;
+
+    }
+
+    if(command_option_exists(argv, argv + argc, "-output_steps"))
     {
-        result.store_delta = true;
+        result.output_steps = true;
     }
 
     return result;

examples/Example_get_apr.h

@@ -6,7 +6,6 @@
 
 #include "algorithm/APRParameters.hpp"
 #include "data_structures/Mesh/PixelData.hpp"
-#include "algorithm/APRConverter.hpp"
 #include "data_structures/APR/APR.hpp"
 
 
@@ -20,8 +19,10 @@ struct cmdLineOptions{
     std::string mask_file = "";
     bool stats_file = false;
     bool normalize_input = false;
-    bool store_delta = false;
+    bool neighborhood_optimization = true;
+    bool output_steps = false;
     unsigned int compress_level = 2;
+    unsigned int compress_type = 0;
 
     float Ip_th = -1;
     float SNR_min = -1;

examples/Example_ray_cast.h

@@ -1,7 +1,6 @@
 #ifndef PARTPLAY_RAYCAST_H
 #define PARTPLAY_RAYCAST_H
 
-
 #include "numerics/APRRaycaster.hpp"
 #include "data_structures/APR/APR.hpp"
 

examples/Example_reconstruct_patch.cpp

@@ -39,6 +39,7 @@ Default: Piece-wise constant reconstruction
 
 #include "data_structures/APR/APR.hpp"
 #include "io/TiffUtils.hpp"
+#include "numerics/APRTreeNumerics.hpp"
 
 
 struct cmdLineOptions{
@@ -186,18 +187,25 @@ int main(int argc, char **argv) {
     reconPatch.z_begin = options.z_begin;
     reconPatch.z_end = options.z_end;
 
+    reconPatch.level_delta = options.level_delta;
+
+    APRTree<uint16_t> aprTree;
+    aprTree.init(apr);
+
     // Intentionaly block-scoped since local recon_pc will be destructed when block ends and release memory.
     {
 
         if(options.output_pc_recon) {
             //create mesh data structure for reconstruction
             PixelData<uint16_t> recon_pc;
 
+            ExtraParticleData<uint16_t> partsTree;
 
+            APRTreeNumerics::fill_tree_from_particles(apr,aprTree,apr.particles_intensities,partsTree,[] (const uint16_t& a,const uint16_t& b) {return std::max(a,b);});
 
             timer.start_timer("pc interp");
             //perform piece-wise constant interpolation
-            aprReconstruction.interp_image_patch(apr,recon_pc, apr.particles_intensities,reconPatch);
+            aprReconstruction.interp_image_patch(apr,aprTree,recon_pc, apr.particles_intensities,partsTree,reconPatch);
             timer.stop_timer();
 
             float elapsed_seconds = timer.t2 - timer.t1;
@@ -210,64 +218,39 @@ int main(int argc, char **argv) {
         }
     }
 
-    //////////////////////////
-    /// Create a particle dataset with the particle type and pc construct it
-    ////////////////////////////
-
-    if(options.output_spatial_properties) {
-
-        //initialization of the iteration structures
-        APRIterator<uint16_t> apr_iterator(apr); //this is required for parallel access
-
-        //create particle dataset
-        ExtraParticleData<uint16_t> type(apr);
-        ExtraParticleData<uint16_t> level(apr);
-
-        ExtraParticleData<uint16_t> x(apr);
-        ExtraParticleData<uint16_t> y(apr);
-        ExtraParticleData<uint16_t> z(apr);
 
-        timer.start_timer("APR parallel iterator loop");
-#ifdef HAVE_OPENMP
-#pragma omp parallel for schedule(static) firstprivate(apr_iterator)
-#endif
-        for (uint64_t particle_number = 0; particle_number < apr_iterator.total_number_particles(); ++particle_number) {
-            //needed step for any parallel loop (update to the next part)
-            apr_iterator.set_iterator_to_particle_by_number(particle_number);
-            type[apr_iterator] = apr_iterator.type();
-            level[apr_iterator] = apr_iterator.level();
+    {
 
-            x[apr_iterator] = apr_iterator.x();
-            y[apr_iterator] = apr_iterator.y();
-            z[apr_iterator] = apr_iterator.z();
-        }
-        timer.stop_timer();
+        if(options.output_smooth_recon) {
+            //create mesh data structure for reconstruction
+            PixelData<uint16_t> recon_pc;
 
-        // Intentionaly block-scoped since local type_recon will be destructed when block ends and release memory.
-        {
-            PixelData<uint16_t> type_recon;
+            ExtraParticleData<uint16_t> partsTree;
 
-            aprReconstruction.interp_image_patch(apr,type_recon, type,reconPatch);
-            TiffUtils::saveMeshAsTiff(options.directory + apr.name + "_type.tif", type_recon);
+            std::vector<float> scale = {1,1,1};
 
-            //level
-            aprReconstruction.interp_image_patch(apr,type_recon, level,reconPatch);
-            TiffUtils::saveMeshAsTiff(options.directory + apr.name + "_level.tif", type_recon);
+            APRTreeNumerics::fill_tree_from_particles(apr,aprTree,apr.particles_intensities,partsTree,[] (const uint16_t& a,const uint16_t& b) {return std::max(a,b);});
 
-            //x
-            aprReconstruction.interp_image_patch(apr,type_recon, x,reconPatch);
-            TiffUtils::saveMeshAsTiff(options.directory + apr.name + "_x.tif", type_recon);
+            timer.start_timer("smooth interp");
+            //perform smooth interpolation
+            aprReconstruction.interp_parts_smooth_patch(apr,aprTree,recon_pc, apr.particles_intensities,partsTree,reconPatch,scale);
+            timer.stop_timer();
 
-            //y
-            aprReconstruction.interp_image_patch(apr,type_recon, y,reconPatch);
-            TiffUtils::saveMeshAsTiff(options.directory + apr.name + "_y.tif", type_recon);
+            float elapsed_seconds = timer.t2 - timer.t1;
+            std::cout << "Smooth recon "
+                      << (recon_pc.x_num * recon_pc.y_num * recon_pc.z_num * 2) / (elapsed_seconds * 1000000.0f)
+                      << " MB per second" << std::endl;
 
-            //z
-            aprReconstruction.interp_image_patch(apr,type_recon, z,reconPatch);
-            TiffUtils::saveMeshAsTiff(options.directory + apr.name + "_z.tif", type_recon);
+            //write output as tiff
+            TiffUtils::saveMeshAsTiff(options.directory + apr.name + "_smooth.tif", recon_pc);
         }
     }
 
+    //////////////////////////
+    /// Create a particle dataset with the particle type and pc construct it
+    ////////////////////////////
+
+
     if(options.output_level) {
 
         //initialization of the iteration structures
@@ -294,7 +277,7 @@ int main(int argc, char **argv) {
             PixelData<uint8_t> type_recon;
 
             //level
-            aprReconstruction.interp_image_patch(apr,type_recon, level,reconPatch);
+            //aprReconstruction.interp_image_patch(apr,aprTree,type_recon, level,reconPatch);
             TiffUtils::saveMeshAsTiff(options.directory + apr.name + "_level.tif", type_recon);
 
         }