get_multiapr can now repeat processing of input data, skip saving APRs and limit log messages for benchmarking purposes

Author: krzysg

examples/Example_get_multiapr.cpp

@@ -12,6 +12,11 @@ Additional settings (High Level):
 -sigma_th   lower threshold for the local intensity scale
 -grad_th    ignore areas in the image where the gradient magnitude is lower than this value
 
+-skipOutputMessages     produce less output/debug messages
+-doNotSaveAPRs          do not save output APR files (good for benchmarking)
+-r                      number of repetitions - default value = 1 means that all input files are processed only once
+                        for higher number input files are processed multiple times (good for benchmarking)
+
 Advanced (Direct) Settings:
 ===========================
 -lambda lambda_value (directly set the value of the gradient smoothing parameter lambda (reasonable range 0.1-10, default: 3)
@@ -43,6 +48,11 @@ struct cmdLineOptions {
     float rel_error = 0.1;
 
     bool neighborhood_optimization = true;
+
+    int numOfRepetitions = 1;
+    int numOfStreams = 3;
+    bool doNotSaveAPRs = false;
+    bool skipOutputMessages = false;
 };
 
 bool command_option_exists(const char **begin, const char **end, const std::string &option)
@@ -112,10 +122,26 @@ cmdLineOptions read_command_line_options(const int argc, const char **argv) {
         options.neighborhood_optimization = false;
     }
 
+    if (command_option_exists(argv, argv + argc, "-skipOutputMessages")) {
+        options.skipOutputMessages = true;
+    }
+
+    if (command_option_exists(argv, argv + argc, "-doNotSaveAPRs")) {
+        options.doNotSaveAPRs = true;
+    }
+
+    if (command_option_exists(argv, argv + argc, "-r")) {
+        options.numOfRepetitions = std::stoi(std::string(get_command_option(argv, argv + argc, "-r")));
+    }
+
+    if (command_option_exists(argv, argv + argc, "-numOfStreams")) {
+        options.numOfStreams = std::stoi(std::string(get_command_option(argv, argv + argc, "-numOfStreams")));
+    }
+
     return options;
 }
 
-
+/* Finds all tiff files (with possible different extensions) in provided directory */
 auto getTiffFilesFromDir(const std::string &directory_path) {
     namespace fs = std::filesystem;
 
@@ -124,8 +150,7 @@ auto getTiffFilesFromDir(const std::string &directory_path) {
     try {
         for (const auto& entry : fs::directory_iterator(directory_path)) {
             if (entry.is_regular_file()) {
-                auto ext = entry.path().extension().string();
-                if (ext == ".tif" || ext == ".tiff" || ext == ".TIF" || ext == ".TIFF") {
+                if (auto ext = entry.path().extension().string(); ext == ".tif" || ext == ".tiff" || ext == ".TIF" || ext == ".TIFF") {
                     tif_files.push_back(entry.path());
                 }
             }
@@ -191,44 +216,61 @@ int runAPR(const cmdLineOptions &options) {
         partIntensities.push_back(std::make_unique<VectorData<ImgType>>(VectorData<ImgType>{}));
         partIntensities_raw.push_back(partIntensities.back().get());
     }
+    // To proces input image multiple times (for benchmarking etc.) we 'multiply input data' by adding extra APR and particle intensity objects and
+    // by copying input raw pointer to images to 'pretend' that we have a lot of input images.
+    size_t numOfInputImages = input_images_raw.size();
+    for (int m = 1; m < options.numOfRepetitions; m++) {
+        for (size_t n = 0; n < numOfInputImages; n++) {
+            input_images_raw.push_back(input_images[n].get());
+            APRs.push_back(std::make_unique<APR>(APR{}));
+            APRs_raw.push_back(APRs.back().get());
+            partIntensities.push_back(std::make_unique<VectorData<ImgType>>(VectorData<ImgType>{}));
+            partIntensities_raw.push_back(partIntensities.back().get());
+        }
+    }
 
     std::cout << std::endl;
 
-    APRTimer timer(true);
+    APRTimer timer(false);
     timer.start_timer("GPU pipeline (mem allocation, processing, sampling) ");
-    if (aprConverter.get_apr_cuda_multistreams(APRs_raw, input_images_raw, partIntensities_raw)) {
+    if (aprConverter.get_apr_cuda_multistreams(APRs_raw, input_images_raw, partIntensities_raw, options.numOfStreams)) {
         timer.stop_timer();
-        size_t numOfImages = input_images_raw.size();
-        std::cout << std::endl;
+        size_t numOfImagesToProcess = input_images_raw.size(); // might be 'processInputMultipleTimes' times bigger than num of input images
+        if (!options.skipOutputMessages) std::cout << std::endl;
 
-        for (size_t i = 0; i < numOfImages; i++) {
-            std::cout << "Postprocessing " << i+1 << "/" << numOfImages << " image...\n";
+        for (size_t i = 0; i < numOfImagesToProcess && !options.doNotSaveAPRs; i++) {
+            if (!options.skipOutputMessages) std::cout << "Postprocessing " << i+1 << "/" << numOfImagesToProcess << " image...\n";
             auto &apr = *APRs[i].get(); // currently process APR
             auto &particle_intensities = *partIntensities[i].get(); // intensities sampled for current APR
 
 
             // ------------ TODO: remove me later, this is quick test for Cpu vs Gpu before real test is written
             // std::cout << apr.linearAccess.y_vec.size() << " particles in APR" << std::endl;
             // std::cout << particle_intensities.size() << " intensities in CPU in APR" << std::endl;
-            // if (apr.linearAccess.y_vec.size() != particle_intensities.size()) {std::cerr << "CPU vs GPU number of particles differ!" << std::endl;}
+            if (apr.linearAccess.y_vec.size() != particle_intensities.size()) {std::cerr << "CPU vs GPU number of particles differ!" << std::endl;}
             ParticleData<ImgType> particle_intensities_cpu;
-            particle_intensities_cpu.sample_image(apr, *input_images[i].get()); // sample your particles from your image
+            particle_intensities_cpu.sample_image(apr, *input_images_raw[i]); // sample your particles from your image
+            int errorCnt = 0;
             for (size_t j = 0 ; j < particle_intensities.size(); ++j) {
                 if (particle_intensities_cpu[j]  != particle_intensities[j]) {
-                    std::cout << "Mismatch at " << j << " CPU: " << particle_intensities_cpu[j] << " GPU: " << particle_intensities[j] << std::endl;
+                    errorCnt++;
+                    // std::cout << "Mismatch at " << j << " CPU: " << particle_intensities_cpu[j] << " GPU: " << particle_intensities[j] << std::endl;
                 }
             }
+            if (errorCnt > 0) std::cout << errorCnt << " errors for index=" << i << std::endl;
             // ---------------------------------------------------------------------------------------------------
 
             // Output name is like base of input filename + extension ".apr"
+            // Extra number is added for 'multiplied' input images
             auto outputDir = std::filesystem::path(options.output_dir);
-            const std::filesystem::path& p(tifFiles[i]);
-            std::string outpuFileName = p.stem().string() + ".apr";
+            const std::filesystem::path& p(tifFiles[i % numOfInputImages]);
+            std::string num = (i >= numOfInputImages) ? std::to_string(i) : "";
+            std::string outputFileName = p.stem().string() + num + ".apr";
 
             //write the APR to hdf5 file
             timer.start_timer("writing output");
             APRFile aprFile;
-            aprFile.open(outputDir / outpuFileName);
+            aprFile.open(outputDir / outputFileName);
             aprFile.write_apr(apr, 0, "t", false);
             ParticleData<ImgType> pd;
             pd.data = std::move(particle_intensities);
@@ -239,10 +281,13 @@ int runAPR(const cmdLineOptions &options) {
             float aprImageSizeInMB = aprFile.current_file_size_MB();
             double originalImageSizeInMB = sizeof(ImgType) * static_cast<double>(apr.org_dims(0) * apr.org_dims(1) * apr.org_dims(2)) / 1'000'000.0;
 
-            std::cout << "Computational Ratio (Pixels/Particles): " << apr.computational_ratio() << std::endl;
-            std::cout << "Original / APR image size:              " << originalImageSizeInMB << " / " << aprImageSizeInMB <<" MB" << std::endl;
-            std::cout << "Lossy Compression Ratio:                " << originalImageSizeInMB/aprImageSizeInMB << std::endl;
-            std::cout << std::endl;
+            if (!options.skipOutputMessages) {
+                std::cout << "Save filename: [" << outputFileName << "]" << std::endl;
+                std::cout << "Computational Ratio (Pixels/Particles): " << apr.computational_ratio() << std::endl;
+                std::cout << "Original / APR image size:              " << originalImageSizeInMB << " / " << aprImageSizeInMB <<" MB" << std::endl;
+                std::cout << "Lossy Compression Ratio:                " << originalImageSizeInMB/aprImageSizeInMB << std::endl;
+                std::cout << std::endl;
+            }
         }
     }
     else {

src/algorithm/APRConverter.hpp

@@ -427,40 +427,39 @@ inline bool APRConverter<ImageType>::get_apr_cuda_multistreams(std::vector<APR*>
     // Reduce number of streams to number of images if there are fewer images than streams
     if (numOfImages < numOfStreams) numOfStreams = numOfImages;
 
-    // Use first image to initialize the APR - all other images should have the same dimensions
-    auto input_image = input_images[0];
-
     // Initialize APRs and memory for the pipeline
     for (auto apr : aAPRs) {
-        if (!initPipelineAPR(*apr, *input_image)) return false;
+        // Use first image to initialize the APR - all other images should have the same dimensions
+        if (!initPipelineAPR(*apr, *input_images[0])) return false;
     }
 
-    // Create a temporary image for each stream
-    std::vector<PixelData<ImageType>> tempImages;
+    // Create a pinned buffer which will be linked to GpuProcessingTask handling each stream
+    // These buffers are used to transmit images from CPU to GPU (first input image need to be copied there)
+    std::vector<PixelData<ImageType>> pinnedBuffers;
     std::cout << "Allocating memory for " << numOfStreams << " streams." << std::endl;
     for (int i = 0; i < numOfStreams; ++i) {
-        tempImages.emplace_back(PixelData<T>(*input_image, true /* copy */, true /* pinned memory */));
+        pinnedBuffers.emplace_back(PixelData<T>(*input_images[i], false /* copy */, true /* pinned memory */));
     }
 
      /////////////////////////////////
     /// Pipeline
     /////////////////////////////////
     APRTimer t(true);
 
-    // Create GpuProcessingTask for each stream
+    // Create GpuProcessingTask for each stream and link it with pinnedBuffer
     std::vector<GpuProcessingTask<ImageType>> gpts;
     t.start_timer("Creating GPTS");
     std::vector<std::future<void>> gpts_futures; gpts_futures.resize(numOfStreams);
     for (int i = 0; i < numOfStreams; ++i) {
-        gpts.emplace_back(GpuProcessingTask<ImageType>(tempImages[i], par, aAPRs[0]->level_max()));
+        gpts.emplace_back(GpuProcessingTask<ImageType>(pinnedBuffers[i], par, aAPRs[0]->level_max()));
     }
     t.stop_timer();
 
-
     t.start_timer("GPU processing...");
     // Saturate all the streams with first images
     for (int i = 0; i < numOfStreams; ++i) {
         std::cout << "Processing image " << i << " on stream " << i  << std::endl;
+        pinnedBuffers[i].copyFromMesh(*input_images[i]);
         gpts_futures[i] = std::async(std::launch::async, &GpuProcessingTask<ImageType>::processOnGpu, &gpts[i]);
     }
 
@@ -476,7 +475,7 @@ inline bool APRConverter<ImageType>::get_apr_cuda_multistreams(std::vector<APR*>
         // We have 'numOfImages - numOfStreams' left to process after saturating the streams with first images
         if (s  < numOfImages - numOfStreams) {
             int imageToProcess = s + numOfStreams;
-            tempImages[streamNum].copyFromMesh(*input_images[imageToProcess]);
+            pinnedBuffers[streamNum].copyFromMesh(*input_images[imageToProcess]);
             std::cout << "Processing image " << imageToProcess << " on stream " << streamNum << std::endl;
             gpts_futures[streamNum] = std::async(std::launch::async, &GpuProcessingTask<ImageType>::processOnGpu, &gpts[streamNum]);
         }
@@ -489,8 +488,8 @@ inline bool APRConverter<ImageType>::get_apr_cuda_multistreams(std::vector<APR*>
         aAPRs[s]->apr_initialized = true;
         if (intensities[s] != nullptr) *intensities[s] = std::move(linearAccessGpu.parts);
     }
-
     auto allT = t.stop_timer();
+
     float tpi = allT / (numOfImages);
     std::cout << "Num of images processed: " << numOfImages << "\n";
     std::cout << "Time per image: " << tpi << " seconds\n";

src/algorithm/ComputeGradientCuda.cu

@@ -401,6 +401,11 @@ template <typename U>
 template <typename ImgType>
 class GpuProcessingTask<U>::GpuProcessingTaskImpl {
 
+
+
+    int cudaDevID;
+    cudaError_t cudaerr;
+
     CudaStream cudaStream;
     const cudaStream_t iStream;
 
@@ -465,6 +470,8 @@ class GpuProcessingTask<U>::GpuProcessingTaskImpl {
 public:
 
     GpuProcessingTaskImpl(const PixelData<ImgType> &inputImage, const APRParameters &parameters, int maxLevel) :
+        cudaDevID(0),  // used for now only in testing new ideas, ideally we should discover CPU/GPU architecture and assing corrct GPU to correct CPU which involves also correct memory allocation and thread affinity....
+        cudaerr(cudaSetDevice(cudaDevID)),
         iCpuImage(inputImage),
         iStream(cudaStream.get()),
         image (inputImage, iStream),
@@ -524,7 +531,7 @@ public:
         // Calculate number of blocks to saturate whole SMs
         // Multiply it by 8 to have more smaller blocks to have better load balancing in case GPU is busy with other tasks
         cudaDeviceProp deviceProp;
-        cudaGetDeviceProperties(&deviceProp, 0);
+        cudaGetDeviceProperties(&deviceProp, cudaDevID);
         const int smCount = deviceProp.multiProcessorCount;
         const int numOfThreadsPerSM = deviceProp.maxThreadsPerMultiProcessor;
         constexpr int numOfThreads = 512;
@@ -563,10 +570,12 @@ public:
     }
 
     LinearAccessCudaStructs<ImgType> getDataFromGpu() {
+        cudaSetDevice(cudaDevID);
         return std::move(lacs);
     }
 
     void processOnGpu() {
+        cudaSetDevice(cudaDevID);
         // Set it and copy first before copying the image
         // It improves *a lot* performance even though it is needed later in computeLinearStructureCuda()
         iAprInfo.total_number_particles = 0; // reset total_number_particles to 0

src/data_structures/Mesh/PixelData.hpp

@@ -19,6 +19,7 @@
 #include <iomanip>
 #include <algorithm>
 #include <numeric>
+#include <cstring>
 
 #include "misc/APRTimer.hpp"
 
@@ -435,8 +436,12 @@ public :
     PinnedMemoryUniquePtr<T> meshMemoryPinned;
 #endif
     ArrayWrapper<T> mesh;
-    
-    uint64_t size() { return (uint64_t) x_num * y_num * z_num * sizeof(T); }
+
+    /**
+     * @return Size in bytes of PixelData (number of elements * sizeof of element)
+     */
+    uint64_t size() const { return (uint64_t) x_num * y_num * z_num * sizeof(T); }
+
     /**
      * Constructor - initialize mesh with size of 0,0,0
      */
@@ -598,7 +603,16 @@ public :
     }
 
     /**
-     * Copies data from aInputMesh utilizing parallel copy, requires prior initialization
+     * Copies data from aInputMesh - requires prior initialization
+     * of 'this' object (size and number of elements)
+     * @param aInputMesh input mesh with data
+    */
+    void copyFromMesh(const PixelData<T> &aInputMesh) {
+        std::memcpy(mesh.begin(), aInputMesh.mesh.begin(), aInputMesh.size());
+    }
+
+    /**
+     * Copies data from aInputMesh which is of different type than 'this' utilizing parallel copy, requires prior initialization
      * of 'this' object (size and number of elements)
      * @tparam U type of data
      * @param aInputMesh input mesh with data

test/MeshDataTest.cpp

@@ -423,6 +423,21 @@ namespace {
         ASSERT_EQ(m(yLen, xLen, zLen), valueForIndex(yLen-1, xLen-1, zLen-1));
     }
 
+    TEST_F(MeshDataTest, CopyFromMeshTest) {
+        PixelData<decltype(m)::value_type> mNew(yLen, xLen, zLen);
+
+        mNew.copyFromMesh(m);
+
+        // Compare if same
+        for (int y = 0; y < yLen; ++y) {
+            for (int x = 0; x < xLen; ++x) {
+                for (int z = 0; z < zLen; ++z) {
+                    // ASSERT_EQ(m(y, x, z), mNew(y, x, z));
+                }
+            }
+        }
+    }
+
     TEST_P(MeshDataParameterTest, BlockCopyDataTest) {
         PixelData<unsigned short> mNew(yLen, xLen, zLen);