Multiple minor fixes and updates: Updates the apr compress method, making go across all particle levels and removing bug. Removed dependency of BLOSC from writing paraview/matlab readable files. Added more file size output to the file IO

cheesema · cheesema · commit e6bbed97b0dd · 2018-05-29T15:42:55.000+02:00
diff --git a/examples/Example_compress_apr.cpp b/examples/Example_compress_apr.cpp
@@ -16,6 +16,7 @@ Example_compress_apr -i input_image_tiff -d input_directory
 
 -compress_type number (1 or 2) (1 - WNL compression (Default), 2 - prediction step with lossless, potential rounding error)
 -quantization_level (Default 1: higher increasing the loss nature of the WNL compression aproach)
+-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.)
 
 e.g. Example_compress_apr -i nuc_apr.h5 -d /Test/Input_examples/ -compress_type 2
 
@@ -59,17 +60,47 @@ int main(int argc, char **argv) {
     APRCompress<uint16_t> comp;
     ExtraParticleData<uint16_t> symbols;
 
+    //feel free to change
+    unsigned int blosc_comp_type = BLOSC_ZSTD;
+    unsigned int blosc_comp_level = options.compress_level;
+    unsigned int blosc_shuffle = 1;
+
     comp.set_quantization_factor(options.quantization_level); //set this to adjust the compression factor for WNL
     comp.set_compression_type(options.compress_type);
 
-    timer.start_timer("compress");
-    apr.write_apr(options.directory ,name + "_compress",comp,BLOSC_ZSTD,1,2);
+    //compress the APR and write to disk
+    timer.start_timer("compress and write");
+    FileSizeInfo fileSizeInfo = apr.write_apr(options.directory ,name + "_compress",comp,blosc_comp_type,blosc_comp_level,blosc_shuffle);
     timer.stop_timer();
 
-    timer.start_timer("decompress");
+    float time_write = (float) timer.timings.back();
+
+    //read the APR and decompress
+    timer.start_timer("read and decompress");
     apr.read_apr(options.directory + name + "_compress_apr.h5");
     timer.stop_timer();
 
+    float time_read = (float) timer.timings.back();
+
+    float original_pixel_image_size = (2.0f*apr.orginal_dimensions(0)*apr.orginal_dimensions(1)*apr.orginal_dimensions(2))/(1000000.0);
+    std::cout << std::endl;
+    std::cout << std::endl;
+    std::cout << "Original image size: " << original_pixel_image_size << " MB" << std::endl;
+
+    float apr_compressed_file_size = fileSizeInfo.total_file_size;
+
+    std::cout << "Compressed (Lossy - WNL) APR: " << apr_compressed_file_size << " MB" << std::endl;
+    std::cout << "Compression Ratio: " << original_pixel_image_size/apr_compressed_file_size << std::endl;
+    std::cout << std::endl;
+    std::cout << std::endl;
+
+    std::cout << "Effective Datarate Write (by original image size): " << original_pixel_image_size/time_write << " MB*/s" << std::endl;
+    std::cout << "Effective Datarate Read (by original image size): " << original_pixel_image_size/time_read << " MB*/s" << std::endl;
+
+    std::cout << std::endl;
+    std::cout << std::endl;
+
+    //writes the piece-wise constant reconstruction of the APR to file for comparison
     PixelData<uint16_t> img;
     apr.interp_img(img,apr.particles_intensities);
     std::string output = options.directory + name + "_compress.tif";
@@ -93,49 +124,53 @@ char* get_command_option(char **begin, char **end, const std::string &option)
 }
 
 
-    cmdLineOptions read_command_line_options(int argc, char **argv){
+cmdLineOptions read_command_line_options(int argc, char **argv){
 
-        cmdLineOptions result;
+    cmdLineOptions result;
 
-        if(argc == 1) {
-            std::cerr << usage << std::endl;
-            exit(1);
-        }
+    if(argc == 1) {
+        std::cerr << usage << std::endl;
+        exit(1);
+    }
 
-        if(command_option_exists(argv, argv + argc, "-i"))
-        {
-            result.input = std::string(get_command_option(argv, argv + argc, "-i"));
-        } else {
-            std::cout << "Input file required" << std::endl;
-            exit(2);
-        }
+    if(command_option_exists(argv, argv + argc, "-i"))
+    {
+        result.input = std::string(get_command_option(argv, argv + argc, "-i"));
+    } else {
+        std::cout << "Input file required" << std::endl;
+        exit(2);
+    }
 
-        if(command_option_exists(argv, argv + argc, "-d"))
-        {
-            result.directory = std::string(get_command_option(argv, argv + argc, "-d"));
-        }
+    if(command_option_exists(argv, argv + argc, "-d"))
+    {
+        result.directory = std::string(get_command_option(argv, argv + argc, "-d"));
+    }
 
-        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, "-compress_type"))
+    {
+        result.compress_type = (unsigned int)std::stoi(std::string(get_command_option(argv, argv + argc, "-compress_type")));
+    }
 
-        if(result.compress_type > 2 || result.compress_type == 0){
+    if(result.compress_type > 2 || result.compress_type == 0){
 
-            std::cerr << "Invalid Compression setting (1 or 2)" << std::endl;
-            exit(1);
-        }
+        std::cerr << "Invalid Compression setting (1 or 2)" << std::endl;
+        exit(1);
+    }
 
 
-        if(command_option_exists(argv, argv + argc, "-quantization_level"))
-        {
-            result.quantization_level =std::stof(std::string(get_command_option(argv, argv + argc, "-quantization_level")));
-        }
+    if(command_option_exists(argv, argv + argc, "-quantization_level"))
+    {
+        result.quantization_level =std::stof(std::string(get_command_option(argv, argv + argc, "-quantization_level")));
+    }
 
+    if(command_option_exists(argv, argv + argc, "-compress_level"))
+    {
+        result.compress_level = (unsigned int)std::stoi(std::string(get_command_option(argv, argv + argc, "-compress_level")));
+    }
 
 
-        return result;
+    return result;
 
-    }
+}
 
 
diff --git a/examples/Example_compress_apr.h b/examples/Example_compress_apr.h
@@ -15,6 +15,7 @@ struct cmdLineOptions{
     unsigned int compress_type=1;
     float quantization_level=1;
     bool stats_file = false;
+    unsigned int compress_level = 2;
 };
 
 cmdLineOptions read_command_line_options(int argc, char **argv);
diff --git a/examples/Example_get_apr.cpp b/examples/Example_get_apr.cpp
@@ -33,8 +33,6 @@ Advanced (Direct) Settings:
 #include "ConfigAPR.h"
 #include "Example_get_apr.h"
 
-
-
 int main(int argc, char **argv) {
 
     //input parsing
diff --git a/src/data_structures/APR/APR.hpp b/src/data_structures/APR/APR.hpp
@@ -55,15 +55,15 @@ class APR {
         apr_writer.read_apr(*this,file_name);
     }
 
-    float write_apr(std::string save_loc,std::string file_name){
+    FileSizeInfo write_apr(std::string save_loc,std::string file_name){
         return apr_writer.write_apr(*this, save_loc,file_name);
     }
 
-    float write_apr(std::string save_loc,std::string file_name,APRCompress<ImageType>& apr_compressor,unsigned int blosc_comp_type,unsigned int blosc_comp_level,unsigned int blosc_shuffle){
+    FileSizeInfo write_apr(std::string save_loc,std::string file_name,APRCompress<ImageType>& apr_compressor,unsigned int blosc_comp_type,unsigned int blosc_comp_level,unsigned int blosc_shuffle){
         return apr_writer.write_apr((*this),save_loc, file_name, apr_compressor,blosc_comp_type ,blosc_comp_level,blosc_shuffle);
     }
 
-    float write_apr(std::string save_loc,std::string file_name,unsigned int blosc_comp_type,unsigned int blosc_comp_level,unsigned int blosc_shuffle){
+    FileSizeInfo write_apr(std::string save_loc,std::string file_name,unsigned int blosc_comp_type,unsigned int blosc_comp_level,unsigned int blosc_shuffle){
         APRCompress<ImageType> apr_compressor;
         apr_compressor.set_compression_type(0);
         return apr_writer.write_apr((*this),save_loc, file_name, apr_compressor,blosc_comp_type ,blosc_comp_level,blosc_shuffle);
@@ -151,9 +151,9 @@ class APR {
         APRIterator<ImageType> apr_iterator(*this); //this is required for parallel access
         parts.data.resize(apr_iterator.total_number_particles());
 
-        #ifdef HAVE_OPENMP
-	    #pragma omp parallel for schedule(static) firstprivate(apr_iterator)
-        #endif
+#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);
diff --git a/src/io/APRWriter.hpp b/src/io/APRWriter.hpp
@@ -13,6 +13,13 @@
 #include <memory>
 
 
+struct FileSizeInfo {
+    float total_file_size=0;
+    float intensity_data=0;
+    float access_data=0;
+};
+
+
 struct AprType {hid_t hdf5type; const char * const typeName;};
 namespace AprTypes  {
 
@@ -174,7 +181,7 @@ class APRWriter {
     }
 
     template<typename ImageType>
-    float write_apr(APR<ImageType>& apr, const std::string &save_loc, const std::string &file_name) {
+    FileSizeInfo write_apr(APR<ImageType>& apr, const std::string &save_loc, const std::string &file_name) {
         APRCompress<ImageType> apr_compressor;
         apr_compressor.set_compression_type(0);
         return write_apr(apr, save_loc, file_name, apr_compressor);
@@ -184,13 +191,14 @@ class APRWriter {
      * Writes the APR to the particle cell structure sparse format, using the p_map for reconstruction
      */
     template<typename ImageType>
-    float write_apr(APR<ImageType> &apr, const std::string &save_loc, const std::string &file_name, APRCompress<ImageType> &apr_compressor, unsigned int blosc_comp_type = BLOSC_ZSTD, unsigned int blosc_comp_level = 2, unsigned int blosc_shuffle=1) {
+    FileSizeInfo write_apr(APR<ImageType> &apr, const std::string &save_loc, const std::string &file_name, APRCompress<ImageType> &apr_compressor, unsigned int blosc_comp_type = BLOSC_ZSTD, unsigned int blosc_comp_level = 2, unsigned int blosc_shuffle=1) {
         APRTimer write_timer;
         write_timer.verbose_flag = false;
 
         std::string hdf5_file_name = save_loc + file_name + "_apr.h5";
         AprFile f{hdf5_file_name, AprFile::Operation::WRITE};
-        if (!f.isOpened()) return 0;
+        FileSizeInfo fileSizeInfo1;
+        if (!f.isOpened()) return fileSizeInfo1;
 
         // ------------- write metadata -------------------------
         writeAttr(AprTypes::NumberOfXType, f.groupId, &apr.apr_access.org_dims[1]);
@@ -225,14 +233,7 @@ class APRWriter {
         writeAttr(AprTypes::NoiseSdEstimateType, f.groupId, &apr.parameters.noise_sd_estimate);
         writeAttr(AprTypes::BackgroundIntensityEstimateType, f.groupId, &apr.parameters.background_intensity_estimate);
 
-        // ------------- write data ----------------------------
-        write_timer.start_timer("intensities");
-        if (compress_type_num > 0){
-            apr_compressor.compress(apr,apr.particles_intensities);
-        }
-        hid_t type = Hdf5Type<ImageType>::type();
-        writeData({type, AprTypes::ParticleIntensitiesType}, f.objectId, apr.particles_intensities.data, blosc_comp_type, blosc_comp_level, blosc_shuffle);
-        write_timer.stop_timer();
+
 
         write_timer.start_timer("access_data");
         MapStorageData map_data;
@@ -262,13 +263,35 @@ class APRWriter {
             writeAttr(AprTypes::NumberOfLevelZType, i, f.groupId, &z_num);
         }
 
+
         // ------------- output the file size -------------------
         hsize_t file_size = f.getFileSize();
+        double sizeMB_access = file_size / 1e6;
+
+        FileSizeInfo fileSizeInfo;
+        fileSizeInfo.access_data = sizeMB_access;
+
+        // ------------- write data ----------------------------
+        write_timer.start_timer("intensities");
+        if (compress_type_num > 0){
+            apr_compressor.compress(apr,apr.particles_intensities);
+        }
+        hid_t type = Hdf5Type<ImageType>::type();
+        writeData({type, AprTypes::ParticleIntensitiesType}, f.objectId, apr.particles_intensities.data, blosc_comp_type, blosc_comp_level, blosc_shuffle);
+        write_timer.stop_timer();
+
+        // ------------- output the file size -------------------
+        file_size = f.getFileSize();
         double sizeMB = file_size / 1e6;
-        std::cout << "HDF5 Filesize: " << sizeMB << " MB\n" << "Writing Complete" << std::endl;
-        return sizeMB;
+
+        fileSizeInfo.total_file_size = sizeMB;
+        fileSizeInfo.intensity_data = fileSizeInfo.total_file_size - fileSizeInfo.access_data;
+
+        std::cout << "HDF5 Total Filesize: " << sizeMB << " MB\n" << "Writing Complete" << std::endl;
+        return fileSizeInfo;
     }
 
+
     template<typename ImageType,typename T>
     void write_apr_paraview(APR<ImageType> &apr, const std::string &save_loc, const std::string &file_name, const ExtraParticleData<T> &parts) {
         std::string hdf5_file_name = save_loc + file_name + "_paraview.h5";
@@ -284,10 +307,7 @@ class APRWriter {
         writeAttr(AprTypes::TotalNumberOfParticlesType, f.groupId, &apr.apr_access.total_number_particles);
 
         // ------------- write data ----------------------------
-        unsigned int blosc_comp_level = 1;
-        unsigned int blosc_shuffle = 2;
-        unsigned int blosc_comp_type = BLOSC_ZSTD;
-        writeData({(Hdf5Type<T>::type()), AprTypes::ParticlePropertyType}, f.objectId, parts.data, blosc_comp_type, blosc_comp_level, blosc_shuffle);
+        writeDataStandard({(Hdf5Type<T>::type()), AprTypes::ParticlePropertyType}, f.objectId, parts.data);
 
         APRIterator<ImageType> apr_iterator(apr);
         std::vector<uint16_t> xv(apr_iterator.total_number_particles());
@@ -296,9 +316,9 @@ class APRWriter {
         std::vector<uint8_t> levelv(apr_iterator.total_number_particles());
         std::vector<uint8_t> typev(apr_iterator.total_number_particles());
 
-        #ifdef HAVE_OPENMP
-	    #pragma omp parallel for schedule(static) firstprivate(apr_iterator)
-        #endif
+#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) {
             apr_iterator.set_iterator_to_particle_by_number(particle_number);
             xv[particle_number] = apr_iterator.x_global();
@@ -307,11 +327,11 @@ class APRWriter {
             levelv[particle_number] = apr_iterator.level();
             typev[particle_number] = apr_iterator.type();
         }
-        writeData(AprTypes::ParaviewXType, f.objectId, xv, blosc_comp_type, blosc_comp_level, blosc_shuffle);
-        writeData(AprTypes::ParaviewYType, f.objectId, yv, blosc_comp_type, blosc_comp_level, blosc_shuffle);
-        writeData(AprTypes::ParaviewZType, f.objectId, zv, blosc_comp_type, blosc_comp_level, blosc_shuffle);
-        writeData(AprTypes::ParaviewLevelType, f.objectId, levelv, blosc_comp_type, blosc_comp_level, blosc_shuffle);
-        writeData(AprTypes::ParaviewTypeType, f.objectId, typev, blosc_comp_type, blosc_comp_level, blosc_shuffle);
+        writeDataStandard(AprTypes::ParaviewXType, f.objectId, xv);
+        writeDataStandard(AprTypes::ParaviewYType, f.objectId, yv);
+        writeDataStandard(AprTypes::ParaviewZType, f.objectId, zv);
+        writeDataStandard(AprTypes::ParaviewLevelType, f.objectId, levelv);
+        writeDataStandard(AprTypes::ParaviewTypeType, f.objectId, typev);
 
         // TODO: This needs to be able extended to handle more general type, currently it is assuming uint16
         write_main_paraview_xdmf_xml(save_loc,hdf5_file_name, file_name,apr_iterator.total_number_particles());
@@ -487,6 +507,13 @@ class APRWriter {
         hdf5_write_data_blosc(aObjectId, aType.hdf5type, aType.typeName, rank, dims, aContainer.data(), blosc_comp_type, blosc_comp_level, blosc_shuffle);
     }
 
+    template<typename T>
+    void writeDataStandard(const AprType &aType, hid_t aObjectId, T aContainer) {
+        hsize_t dims[] = {aContainer.size()};
+        const hsize_t rank = 1;
+        hdf5_write_data_standard(aObjectId, aType.hdf5type, aType.typeName, rank, dims, aContainer.data());
+    }
+
     void writeString(AprType aTypeName, hid_t aGroupId, const std::string &aValue) {
         if (aValue.size() > 0){
             hid_t aid = H5Screate(H5S_SCALAR);
diff --git a/src/io/hdf5functions_blosc.cpp b/src/io/hdf5functions_blosc.cpp
@@ -69,6 +69,37 @@ void hdf5_write_data_blosc(hid_t obj_id, hid_t type_id, const char *ds_name, hsi
     H5Pclose(plist_id);
 }
 
+/**
+ * writes data to the hdf5 file or group identified by obj_id of hdf5 datatype data_type without using blosc
+ */
+void hdf5_write_data_standard(hid_t obj_id, hid_t type_id, const char *ds_name, hsize_t rank, hsize_t *dims, void *data) {
+    hid_t plist_id  = H5Pcreate(H5P_DATASET_CREATE);
+
+    // Dataset must be chunked for compression
+    const uint64_t max_size = 100000;
+    hsize_t cdims = (dims[0] < max_size) ? dims[0] : max_size;
+    rank = 1;
+    H5Pset_chunk(plist_id, rank, &cdims);
+
+    //compression parameters
+    int deflate_level = 9;
+
+    /////SET COMPRESSION TYPE /////
+
+    //DEFLATE ENCODING (GZIP)
+    H5Pset_deflate (plist_id, deflate_level);
+
+    //create write and close
+    hid_t space_id = H5Screate_simple(rank, dims, NULL);
+    hid_t dset_id = H5Dcreate2(obj_id, ds_name, type_id, space_id, H5P_DEFAULT, plist_id, H5P_DEFAULT);
+    H5Dwrite(dset_id,type_id,H5S_ALL,H5S_ALL,H5P_DEFAULT,data);
+    H5Dclose(dset_id);
+
+    H5Pclose(plist_id);
+
+}
+
+
 /**
  * writes data to the hdf5 file or group identified by obj_id of hdf5 datatype data_type
  */
diff --git a/src/io/hdf5functions_blosc.h b/src/io/hdf5functions_blosc.h
@@ -28,6 +28,7 @@ void hdf5_load_data_blosc(hid_t obj_id, void* buff, const char* data_name);
 void hdf5_load_data_blosc(hid_t obj_id, hid_t dataType, void* buff, const char* data_name);
 void hdf5_write_attribute_blosc(hid_t obj_id,hid_t type_id,const char* attr_name,hsize_t rank,hsize_t* dims, const void * const data );
 void hdf5_write_data_blosc(hid_t obj_id,hid_t type_id,const char* ds_name,hsize_t rank,hsize_t* dims, void* data ,unsigned int comp_type,unsigned int comp_level,unsigned int shuffle);
+void hdf5_write_data_standard(hid_t obj_id,hid_t type_id,const char* ds_name,hsize_t rank,hsize_t* dims, void* data );
 void write_main_paraview_xdmf_xml(const std::string &aDestinationDir,const std::string &aHdf5FileName, const std::string &aParaviewFileName, uint64_t aNumOfParticles);
 
 
diff --git a/src/numerics/APRCompress.hpp b/src/numerics/APRCompress.hpp
