Merge pull request #91 from cheesema/develop

Develop

Author: cheesema

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_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);
 
         }

src/algorithm/APRConverter.hpp

@@ -27,6 +27,8 @@ template<typename ImageType>
 class APRConverter: public LocalIntensityScale, public ComputeGradient, public LocalParticleCellSet, public PullingScheme {
 
 public:
+
+
     APRParameters par;
     APRTimer fine_grained_timer;
     APRTimer method_timer;
@@ -52,20 +54,21 @@ class APRConverter: public LocalIntensityScale, public ComputeGradient, public L
         }
     };
 
-private:
     //get apr without setting parameters, and with an already loaded image.
     template<typename T>
     bool get_apr_method(APR<ImageType> &aAPR, PixelData<T> &input_image);
 
+    template<typename T>
+    void auto_parameters(const PixelData<T> &input_img);
+
+private:
+
     //pointer to the APR structure so member functions can have access if they need
     const APR<ImageType> *apr;
 
     template<typename T>
     void init_apr(APR<ImageType>& aAPR, PixelData<T>& input_image);
 
-    template<typename T>
-    void auto_parameters(const PixelData<T> &input_img);
-
     template<typename T>
     bool get_apr_method_from_file(APR<ImageType> &aAPR, const TiffUtils::TiffInfo &aTiffFile);
 
@@ -191,15 +194,23 @@ bool APRConverter<ImageType>::get_apr_method(APR<ImageType> &aAPR, PixelData<T>&
     fine_grained_timer.stop_timer();
 
 #ifndef APR_USE_CUDA
-    method_timer.verbose_flag = true;
+    //method_timer.verbose_flag = true;
     method_timer.start_timer("compute_gradient_magnitude_using_bsplines");
     get_gradient(image_temp, grad_temp, local_scale_temp, local_scale_temp2, bspline_offset, par);
     method_timer.stop_timer();
 
+    if(par.output_steps){
+        TiffUtils::saveMeshAsTiff(par.output_dir + "gradient_step.tif", grad_temp);
+    }
+
     method_timer.start_timer("compute_local_intensity_scale");
     get_local_intensity_scale(local_scale_temp, local_scale_temp2, par);
     method_timer.stop_timer();
-    method_timer.verbose_flag = false;
+    //method_timer.verbose_flag = false;
+
+    if(par.output_steps){
+        TiffUtils::saveMeshAsTiff(par.output_dir + "local_intensity_scale_step.tif", local_scale_temp);
+    }
 #else
     method_timer.start_timer("compute_gradient_magnitude_using_bsplines and local instensity scale CUDA");
     getFullPipeline(image_temp, grad_temp, local_scale_temp, local_scale_temp2,bspline_offset, par);
@@ -234,8 +245,6 @@ bool APRConverter<ImageType>::get_apr_method(APR<ImageType> &aAPR, PixelData<T>&
 
     computation_timer.stop_timer();
 
-    aAPR.parameters = par;
-
     total_timer.stop_timer();
 
     return true;
@@ -268,6 +277,12 @@ void APRConverter<ImageType>::get_local_particle_cell_set(PixelData<ImageType> &
     fine_grained_timer.start_timer("compute_level_second");
     //incorporate other factors and compute the level of the Particle Cell, effectively construct LPC L_n
     compute_level_for_array(local_scale_temp,level_factor,par.rel_error);
+
+    if(par.output_steps){
+        TiffUtils::saveMeshAsTiff(par.output_dir + "local_particle_set_level_step.tif", local_scale_temp);
+    }
+
+
     fill(l_max,local_scale_temp);
     fine_grained_timer.stop_timer();
 
@@ -293,7 +308,7 @@ void APRConverter<ImageType>::get_gradient(PixelData<ImageType> &image_temp, Pix
     //  Input: full sized image.
     //  Output: down-sampled by 2 gradient magnitude (Note, the gradient is calculated at pixel level then maximum down sampled within the loops below)
 
-    fine_grained_timer.verbose_flag = true;
+    //fine_grained_timer.verbose_flag = true;
 
     fine_grained_timer.start_timer("threshold");
 
@@ -416,6 +431,8 @@ void APRConverter<ImageType>::init_apr(APR<ImageType>& aAPR,PixelData<T>& input_
 
     aAPR.apr_access.level_min = levelMin;
     aAPR.apr_access.level_max = levelMax;
+
+    aAPR.parameters = par;
 }
 
 template<typename ImageType> template<typename T>